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ABSTRACT 


Yerkes,  Kirk  Lee,  M.S.,  Department  of  Mechanical  and  Materials  Engineer¬ 
ing,  Wright  State  University.  1989.  A  Study  of  Mixed  Convection  in  Large  Baffled 
Rectangular  Chambers  with  and  without  Internal  Heat  Sources. 


A  numerical  and  experimental  investigation  to  determine  the  thermal  effects 
on  the  development  of  the  flow  structure  in  large  baffled  rectangular  chambers 
with  and  without  internal  heat  sources  has  been  completed.  Two-  and  three- 
dimensional  numerical  models  were  formulated  using  the  time-dependent  laminar 
Navier  Stokes  equations  assuming  a  Boussinesq  fluid  with  a  Prandtl  number  of 
0.7.  Experiments  were  conducted  using  a  scaled  down  laboratory  model  and  t lie 
larger,  full  size  chamber.  The  scaled-down  model  simulated  the  full  size  chamber 
enabling  experimental  data  to  be  obtained  and  subsequently  compared  with  nu¬ 
merical  results.  Flow  visualization  experiments  were  conducted  using  the  larger, 
full  si/e  (hamuef  internal  horizontal  baffles  and  heat  sources  were  located  sym¬ 
metrically  about  the  vertical  centerline  axis.  Mixed  convection  without  internal 
heat  source*  for  both  "aided"  and  “opposed”  buoyant  force*  showed  that  the 
development  of  the  flow  structure  wa*  sensitive  to  small  variation*  in  the  tern 
perafure  difference  bet  wren  the  inlet  and  the  vertical  walls.  Mixed  convection 


tit 


with  interna]  heat  sources  showed  the  flow  structure  to  develop  through  a  se¬ 
ries  of  bifurcations  from  steady  state,  to  periodic,  aperiodic  and  Anally  chaotic 
v.-'h  increasing  heat  source  temperature.  Use  of  the  scaled  down  laboratory 
«  x,.t  rimer. t a.  m<  -  tel  as  a:.  m«i.cat..r  for  the  flow  development  m  trie  larger  full- 
si /e  chamber  showed  significant  three  dimensional  effects.  Flow  visualization  in 
the  larger  full-size  chamber  using  Silahydrocarbon  aerosol  droplets  showed  good 
agreement  with  the  two-dimensional  numerical  results. 
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T„  end  wall  temperature,  °C  (see  Fig.  1  of  chapter  1) 

u  dimensional  x-di:ection  velocity,  m/s 

U  non-dimensional  x-direction  velocity,  v/Win 

t;  dimensional  y-direction  velocity,  m/s 

V  non-dimensional  y-diiection  velocity,  chapter  1,  m/s 

V  fluid  velocity  or  droplet  terminal  settling  velocity,  m/s 

V  velocity  vector,  m/s 

VcK  chamber  volume,  ms 

Kno«  magnitude  of  maximum  velocity,  m/s 

Vth  thermophoretic  velocity  vector,  m/s 

V  volume,  m1 

w  dimensional  z-direction  velocity,  m/s 

W  non-dimensional  i-direction  velocity,  w/Wjs 

Win  inlet  vertical  velocity,  m/s 

x  dimensional  horisontal  coordinate,  m 

X  non-dimensional  horizontal  coordinate,  zj D 
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y  dimensional  horizontal  coordinate,  rotameter  ball  position,  m 
Y  non-dimensional  horizontal  coordinate,  y  j L 

Y’  experimental  horizontal  coordinate  axis 

Z  non-dimensional  vertical  coordinate,  zj L 

Greek  letter  symbols 
a  thermal  diffusivity,  m2/s 

3  coefficient  of  thermal  expansion,  3  —  ~ 2  (Jf)i  1/°K 
S  boundary  layer  thickness,  m 

SP  (P2-P1),  N/m2 

ST  {Tis  -  T„),  °C 

e  fractional  change  in  dependent  variable 

$  non-dimensional  time,  tQiff /Vch 

9o  {'I\0 1  -  Tcou)/Tcou 

A  mean  free  path,  m 

p  dynamic  viscosity  or  reference  at  inlet,  kg/m-s 
v  kinematic  viscosity  or  reference  at  inlet,  pjp,  m2/s 
p  fluid  density  or  reference  at  inlet,  kg/m3 
r  time  constant,  s 

r  stress  tensor,  fiW  4-  /i(VY)7,  O7”  ^transpose 

4>  dependent  variable 

dependent  variable  at  specific  sweep  or  time  step 
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VT  temperature  gradient,  °K/m 

Subscripts 

0,1,2  reference  location  for  rotameter  pressure  correction 
o  reference  value  for  Boussinesq  approximation 

atm  atmosphere  reference 

b  rotameter  ball 

d  droplet 

/  fluid 

wire  anemometer  wire  or  film 
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STUDY  INTRODUCTION 


Large,  low  velocity  chambers  {^sl  m  diameter)  have  traditionally  beeD  used 
for  studies  investigating  the  tome  effects  of  various  airborne  compounds.  Use 
of  these  chambers  in  toxicological  studies  have  precluded  detailed  investigations 
utilizing  computational  fluid  dynamics  in  conjunction  with  experimental  data 
in  determining  the  physics  of  chamber  operation  and  performance  in  the  devel¬ 
opment  of  the  chamber  flow  structure.  The  series  of  studies  contained  herein 
address  these  issues  by  utilizing  computational  fluid  dynamics  and  experimental 
data  from  both  a  scaled-down  laboratory  model  and  a  full-size  chamber  typically 
used  for  toxicological  studies.  Internal  heat  sources,  simulating  biological  heat 
sources  have  also  been  included  and  the  effects  on  the  developing  flow  structure 
with  and  without  these  internal  heat  sources  invest:*ated. 

Chapter  1  addresses  the  development  of  the  chamber  flow  structure  with¬ 
out  inteiual  heat  sources  for  a  geometry  typically  used  for  tox;cological  stud¬ 
ies.  Chapter  2  includes  the  addition  of  internal  heat  sources  and  addresses  the 
time-dependent  effects  on  the  development  of  the  chamber  flow  structure.  Both 
chapters  1  and  2  utilize  computational  fluid  dynamics  and  **rperimental  results 
from  the  sealed  down  laboratory  model.  Finally,  chapter  3  utilizes  both  com¬ 
putational  fluid  dynamics  ana  flow  visualization  experiments  conducted  in  the 
full-size  chamber  with  and  without  interna]  heat  sources. 
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Chapter  1 


EXPERIMENTAL  AND  NUMERICAL  SIMULATION  OF  MIXED 
CONVECTION  IN  LARGE  BAFFLED  RECTANGULAR  CHAMBERS 


1.1  Summary 


An  investigation  to  determine  the  effects  of  mixed  convection  on  the  flow 
structure  of  large  baffled  chambers  for  both  "aided”  and  "opposed”  buoyant 
forces  has  been  completed.  Experimental  results  were  reported  for  momentum 
boundary  layer  thickness  and  temperature  distribution  within  the  chamber.  Nu¬ 
merical  results  and  experimental  data  were  obtained  for  comparison.  Two-  and 
three-dimensional  transient  numerical  models  were  formulated  using  a  fully  im¬ 
plicit  finite- difference  scheme  of  the  laminar  NavierStokes  equations.  The  inlet 
and  wall  temperatures  were  held  constant  while  maintaining  a  uniform  inlet 
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velocity.  Two  tiers  of  fom  (two  per  tier)  horizontal  baffles  were  located  Ajrisym- 
metrically  about  the  vertical  axis  with  the  inlet  at  the  top  and  exhaust  at  the 
base  of  the  chamber.  The  difference  between  the  inlet  and  wall  temperatures 
ranged  between  -1.0CC  to  3.3°C  with  a  Prandtl  number,  Pr  =  0.7,  inlet  Reynolds 
number,  Re  -  32  and  235  and  inlet  Rayleigh  number  Ra  --0  -  6.3  x  107  .  It 
was  concluded  that  the  flow  patterns  within  these  chambers  are  predominantly 
buoyant  in  nature  with  asymmetric  behavior.  The  development  of  the  flow  struc¬ 
ture  was  found  to  be  sensitive  to  small  variations  in  the  tempe/aiure  difference 
between  the  inlet  and  wall. 

1.2  Introduction 


Large  baffled  rectangular  chambers  with  low  inlet  velocities  have  traditionally 
been  used  to  study  a  variety  of  chemical  atmospheres,  generated  for  health  effects 
studies,  chemical  reaction  and  chemical  species  formation  studies,  and  aerosol 
characterization  studies.  Fig  1.1a.  Generally,  they  have  a  cross-sectional  diameter 
greater  than  0.5  meter  with  an  inlet  velocity  set  to  maintain  10-15  chamber 
volume  changes  per  hour.  The  inlet  and  exhaust  mt)  be  mounted  vertically  or 
horizontally  with  an  array  of  internal  baffles  included  for  atmosphere  dispersion 
purposes. 

There  have  been  many  investigations  concerning  the  Rayleigh- Benaxd  flow 
in  confined  and  partially  open  enclosures  with  and  without  baffles  or  dividers 
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(Ostrach  ar.d  Austin,  1983;  Bajorek  and  Lloyd,  1982;  Nansteel  and  Greif,  1981; 
Markatos  and  Maiin,  1982).  Such  investigations  have  made  use  of  both  exper¬ 
imental  and  numerical  results  for  information  with  regard  to  the  temperature 
profile  and  flow  structure  in  various  closed  cavity  geometries  with  differing  mag¬ 
nitudes  of  Grashot.  Prandtl,  and  Rayleigh  numbers.  There  have  a'.ro  been  many 
investigations  concerning  mixed  convective  flow  in  pipe  (Tanaka  et  al.,  1987), 
inclined  tubes  ( Choucthury  and  Patankar,  1SSS),  vertical  channels  (Yao,  1983), 
spherical  annuli  (R&madhyar.i  et  al.,  1934),  and  vertical  and  horizontal  plates 
(Rrmachandran  et  al.,  1938;  Osborne  and  Incropera,  1985;  Chiu  and  Rosen- 
berger,  1987;  Chiu  et  al.  1937),  but  there  is  little  information  with  regard  to 
mixed  con  recti  ve  flow  in  large  baffled  chambers. 


Mixed  convection  in  heated  vertical  channels  have  shown  some  indication  of 
the  possible  effects  of  both  ” aided”  and  “opposed”  buoyant  forces  in  larger  ge¬ 
ometries.  Yao  (1933)  demonstrated  that  the  "aided”  buoyant  forces  for  upward 
vertical  channel  flow  tended  to  accelerate  the  boundary  layer  such  that  incoming 
flow  was  drawn  toward  the  channel  wall.  Tanaka  et  al.  (1937)  also  showed  an 
increase  in  the  wall  boundary  lay  velocity  with  a  decrease  in  Reynolds  num¬ 
ber  for  the  case  <>■  “aided”  buoyant  forces  in  a  uniformly  heated  pipe.  These 
investigations  indicate  a  shift  in  the  major  velocity  component  toward  the  pipe 
or  channel  waits  as  buoyant  forces  become  significant.  While  these  investigations 
indicate  some  pose:  ole  clfcct*  of  buoyant  forces,  the  chambers  of  interest  in  ad¬ 
dition  to  in'luding  baffles,  also  have  inlet  Reynolds  numbers  that  are  one  to  two 


*  -v  V*  •f».3L\T.  t.  ■  •  ,  ■  >m- 
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orders  of  magnitude  lower,  resulting  in  a  flow  structure  that  is  predominantly 
buoyant  in  nature. 

Initial  investigations  of  these  chambers  incorporated  flow  visualization  tech¬ 
niques  using  smoke  and  dyes  to  tiriu/e  the  chamber  geometry  so  that  well- 
dispersed  atmospheres  could  be  obtained  (Carpenter  and  Beethe,  1981;  Moss, 
1981).  Quantitative  approaches  to  determine  thi  mean  residence  time  of  a  cham¬ 
ber,  the  rate  of  dispersion  within  a  chamber  of  a  tracer  gas,  and  the  distribution 
within  a  chamber  using  vapors,  droplets,  and  solid  particles  have  also  been  done 
(Hemcrway  et  a].,  1982;  Beethe  et  al..  1979;  Moss,  1982;  Yeh  et  al.,  1986).  These 
investigations,  however,  have  not  addressed  the  effects  of  mixed  convection  on 
the  transport  phenomenon  of  heat  and  mas.  transfer  within  these  chambers.  In¬ 
formation  into  the  diffusion  and  convection  of  mass,  enthalpy,  and  momentum 
within  these  chandlers  is  of  interest  to  better  understand  their  operation  and 
provide  a  basis  for  improvement,  modification  and  future  development. 

A  typical  chamber,  that  may  be  used  for  the  purposes  described,  consists 
of  a  conical  ir.let  intended  to  uniformly  distribute  the  atmosphere  of  interest 
across  the  inlet  plane.  The  exhaust  consists  of  a  manifold  to  achieve  uniform 
flow  characteristics  across  the  exhaust  plane.  Inlet  and  exhaust  flow  rates  are 
controlled  independently  to  maintain  a  slight  negative  chamber  pressure  with 
respect  to  the  ambient  prr»»ure  (2.5  5.1  cm  HjOj. 

This  investigation  considers  a  simplified  version  of  the  aforementioned  cham¬ 
ber.  Of  interest  is  the  development  of  the  flow  structure  due  to  a  temperature 
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gradient  between  the  inlet  and  wall  temperature.  Table  1.1  shows  the  scope  of 
this  investigation  with  both  two-  and  three-dimensional  numerical  formulations 
(cases  #1  -  #6)  and  experimental  results  from  a  scaled-down  laboratory  experi¬ 
mental  model  (case  #7)  intended  to  be  used  as  an  indicator  for  use  in  the  analy¬ 
sis  of  the  larger,  full-size  chamber.  For  this  investigation,  the  three-dimensional 
effects  of  the  experimental  model  were  considered  and  the  corresponding  three- 
dimensional  numerical  solutions  were  used  to  account  for  discrepancies  between 
the  experimental  and  two-dimensional  numerical  results.  The  two-dimensional 
numerical  formulation  is  intended  to  correspond  to  the  larger,  full-size  chamber 
with  a  depth  of  one  meter  and  was  compared  to  the  experimental  results.  The 
vertical  wall  temperature,  Tw,  was  assumed  to  be  constant  due  to  the  thin  wall 
construction.  A  typical  range  of  Reynolds,  Grashof,  Rayleigh  and  Prandtl  num¬ 
bers  at  the  inlet  plane,  assuming  a  uniform  velocity  distribution,  are  also  shown 
in  Table  1.1  where  ST  is  taken  to  be  the  difference  between  the  inlet  tempera¬ 
ture,  Tjs ,  and  the  vertical  wall  temperature,  Tw  ,  with  the  characteristic  length 
chosen  to  be  the  hydraulic  diameter  of  the  chamber  body. 

Initial  verification  of  the  magnitude  of  the  buoyant  forces  compared  to  inertia 
forces  was  determined  by  finding  the  ratio  of  the  buoyant  forces  to  inertia  forces. 
This  was  done  by  evaluating  Gr/Re3  at  the  inlet  plane.  The  buoyant  forces  were 
found  to  dominate  the  inertia  forces  by  as  much  as  three  orders  of  magnitude  with 
a  temperature  difference  of  ST  —  1.0°C.  The  dominate  buoyant  forces,  therefore, 
must  be  considered  as  the  primary  force  in  development  of  the  flow  structure 
within  these  chambers. 
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1.3  Analysis 


Both  two-  and  three-dimensional  numerical  solutions  of  the  incompressible 
time-dependent  laminar  Navier-Stokes  equations  were  considered.  The  time- 
dependent  approach  was  chosen  to  evaluate  the  presence  of  &  steady  periodic 
oscillation.  The  two-  and  three-dimensional  problems  shown  in  Fig.  1.1  were 
considered  with  the  three-dimensional  formulation  simulating  the  experimental 
model.  The  additional  third-  dimensional  wall  is  termed  "side  wall”  while  the 
constamt  temperature  two-  and  three-dimensional  vertical  wall  is  termed  "end 
wall”. 

For  this  investigation,  only  the  dimensional  form  of  the  governing  equations 
were  solved,  thereby  circumventing  the  immediate  need  for  scaling  considerations 
of  the  non-dimensional  parameters  as  is  typically  required  for  buoyant  flow  con¬ 
ditions  (Ostrach  and  Austin,  1988).  The  non-dimensional  form  of  the  governing 
equations  were  derived  using  the  scale  length  and  characteristic  velocity  gener¬ 
ally  chosen  for  mixed  convection  pipe  flow  and  the  resultant  non-dimensional 
parameters  were  compared  using  subsequent  dimensional  solutions. 

The  Boussinesq  approximation  was  used  in  the  vertical  direction  to  account 
for  buoyancy  effects.  Using  this  approximation  the  vertical  body  force  becomes: 


B*  =  PsU  ~  ~  Tin)] 
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where 


for  an  ideal  gas. 


0  =  tTs 


For  this  investigation,  the  reference  values  for  all  of  the  gas  properties  were 
taken  to  be  those  at  the  inlet  plane. 

The  chamber  inlet  cone  was  assumed  to  provide  a  uniform  velocity  distri¬ 
bution  across  the  inlet  plane  with  the  exhaust  manifold  providing  a  velocity 
distribution  resulting  from  a  constant  exhaust  pressure  across  the  outlet  plane. 
Viscous  dissipation  and  pressure  work  were  assumed  to  be  negligible  due  to  the 
low  buoyancy-induced  velocities.  With  these  assumptions,  the  three-dimensional 
governing  equations  become: 


Conservation  of  mass: 


du  dv  diu 

-  A. - 1-  — -  =  0 

dx  dy  dz 


Conservation  of  momentum: 


(  du  du  du  du\  dp+  ( 

(«+“ai  +  vaJ  +  “‘aij“-ar  +  '‘i 


d3  u  d3u  93u> 
dx 2  dy 3  dz3 , 


( dv  dv  dv  dv\  dp+  ( 

'{r,+'‘T,  +  vei  +  wo;)  =  -!7  +  ,‘{ 


d3V  d3V  d2V ' 

dx 3  dy 3  dz3  4 
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dw  dw  dw  dw\  dp+  (  d2w  d2w  d2w 


^U  +  u^+u%+u;a7j  ="a7  +  Ml^+%I  +  ^ 


+Pff^(T  -  TIN) 


Conservation  of  energy: 


far  dT  dT  dT\  k  (d3T  ,  d2T  ,  a3T\ 

(at +“^  +  ’’^  +  "’a7j  =  ^U?+ W  + 


Defining  the  non-dimensional  parameters  as  follows; 


9  =  tQ1N/Vch  =  iWjs/H ,  X  =  x/D,  Y  =  y/L,  Z  =  z/ff, 


t/  =  n/Wiff,  V  =  v/W/n,  W  —  w/WjSy  T  —{Tjh  —  T)/(T//v  —  Tw)t 

P={?+-  K)  / pWIn 

The  non-dimensional  form  of  the  governing  equations  can  be  derived  in  terms 
of  the  aspect  ratios,  A j,  Aj,  and  /lj. 


Conservation  of  mass: 


ou  dv  dw 

A'dX  +  dY*  1  dZ 


(1.6) 
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Conservation  of  momentum: 

dd  \Ai  ax  Ai  dY  dz) 

AjdP  yij  //A3\2d2U  ,  1  d2U  ,  d2u\ 
AidX  +  Re  VVaJ  dX2  +  A\  dY2  ^  dZ2 ) 


dV 

dX 


+ 


1 

At 


+  W 


i  dp  Aif/A2\2d2v  i  d2v  a2v\ 

Ax  dY  +  Re  [\Ai )  3X2  +  A]  dY2  +  dZ2  ) 


(1.8) 


aw  (A3trd\v  i  aw 

39  V^i  dX  Ax  dY 


dP  Ax((A2\2d2W  1  32W  S2W\  GrT’ 

dZ  +  Re  \\Ax)  dX2  +  A]  dY2  +  dZ2  )  +  A3Re2 


(1.9) 


Conservation  of  energy: 


A  i  dX 


1  3T 
+  —  V  — 
Ai  dY 


'  ,dT'\ 

~  +  vv  Jz ) 


Ax  f(A3\2d2T'  1  d2T'  d2T'\ 

Pr  •  Re  [\aJ  3X 2  +  A]  dY2  +  3Z2  ) 


(1.10) 


The  boundary  conditions  were  specified  with  the  assumption  that  there  is  a 
uniform  inlet  velocity  profile  and  the  vertical  end  wall  temperature  is  constant. 
The  side  walls  were  either  at  a  fixed  temperature,  (T‘)  =  1/6,  or  adiabatic, 

(dT' I3x}  -  0,  for  the  three-dimensional  formulation.  Only  the  magnitude  of 
the  inlet  velocity  is  specified.  The  magnitude  of  the  outlet  velocity  is  such  as 
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to  satisfy  the  conservation  of  mass  with  a  constant  reference  pressure  across  the 
outlet  plane.  With  these  considerations  the  boundary  conditions  become: 


X  =  0,0  <  Y  >  1,0  <  Z  >  1,  {dT'jdx)  =  0  or  T'  =  1/6, U  =  0, 

V  =  0,  W  =  0 

X  =  1,0  <  Y  >  1,0  <  Z  >  1,  ( dT'/dx )  =  0  or  T'  =  1/6,  U  =  0, 

V  =  0,  W  =  0 

y  =  o,o  <  x  >  i,o  <  z  >  i.T’  =  \,u  =  o,y  =  o,iv  =  o 

y  =  i,o  <  x  >  i,o  <  z  >  \,t'  =  i,i/  =  o,y  =  o,»v  =  o 

Z  =  0,0  <  X  >  1,0  <  Y  >  1,T‘  =  0, 1/  =  0,V  =  0,VV  =  1 

z  =  i,o<A'>i,o<y>i,c/  =  o,v'  =  o,p  =  o 


The  outlet  boundary  condition  for  temperature  was  considered  to  be  locally 
parabolic  such  that  the  Peclet  number  is  sufficiently  large  so  as  to  exhibit  lo¬ 
cal  one-way  behavior  in  the  axial  direction.  The  baffles  were  treated  as  being 
infinitesimally  thick  with  a  no-slip  boundary  condition.  The  baffle  temperature 
was  specified  to  be  that  of  the  local  fluid  temperature. 
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J.4  Numerical  Scheme 


The  numerical  scheme  used  is  a  finite- difference  iterative  method  of  solution 
using  a  control-volume  approach  as  developed  by  Spalding  et  al.  (1980).  Two- 
and  three-dimensional  numerical  models  were  formulated  using  the  fully  implicit 
finite- difference  scheme  of  the  time  dependent,  laminar  Navier-Stokes  equations. 
The  "SIMPLEST”  (Spalding  et  al.,  1980)  method  of  solution  for  the  momentum 
equations  was  used  with  the  hybrid  differencing  formulation.  When  a  cell  Peclet 
number  is  within  the  range  -2  to  2,  a  central-difference  scheme  is  used  and  when 
the  cell  Peclet  number  is  outside  this  range  the  upwind  differencing  scheme  is 
used. 

The  outlet  boundary  condition  for  temperature  was  considered  to  be  locally 
parabolic  such  that  the  Peclet  number  is  sufficiently  large  so  as  to  exhibit  local 
one-way  behavior  in  the  axial  direction  (Patankar,  1980).  This  approach  to 
specifying  the  outlet  boundary  condition  for  a  dependent  variable  is  generally 
accepted  as  the  most  available  approach  for  internal  flow  problems  in  which 
the  fluid  leaves  the  calculation  domain  and  the  dependent  variable  of  interest  is 
unknown. 

The  baffles  were  treated  as  being  infinitesimally  thick  with  a  no-slip  bound¬ 
ary  condition  by  setting  the  cell  wall  velocity  to  zero  over  the  region  of  the 
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baffle  location.  The  baffle  temperature  was  specified  to  be  that  of  the  local  fluid 
temperature. 

The  solution  sequence  involved  solving  for  the  velocity  and  temperature  fields 
from  an  assumed  pressure  field  at  a  specific  time  step.  The  pre;  .  ire  field  was 
then  subsequently  updated  using  these  velocity  and  temperature  fields  in  the 
pressure-correction  equation  such  that  the  conservation  of  mass  was  satisfied. 
This  iterative  sweep  process  was  then  repeated  using  a  slab-by-slab  method  at 
each  time  step  until  convergence  criteria  were  satisfied.  Convergence  charac¬ 
teristics  at  each  time  step  were  maintained  by  reducing  the  time  step  where 
appropriate  with  the  minimum  computational  time  required  to  reach  a  steady 
state  solution  being  at  least  three  characteristic  times,  /f/W/jv,  or  approximately 
ten  minutes.  In  addition,  conditions  for  convergence  to  a  solution,  either  at  a 
specific  time  step  or  to  a  steady  state  solution,  was  based  upon  a  dependent 
variable,  <j>,  varying  less  than  a  predetermined  change  of  the  magnitude  fraction, 
c^,  between  successive  sweeps  or  time  steps  where: 

I  <£»-<£*  •*•  l  I  /4>  <  (4 

The  values  of  were  varied  and  it  was  required  that  the  magnitude  of  the 
dependent  variable  be  greater  than  the  roundoff  error  or  truncation  error  such 
as  would  happen  when  velocities  approach  zero.  Typically,  for  a  steady  state 
solution,  ranged  from  1.0  x  10" 1  to  1.0  x  10“3  over  a  time  span  of  25  seconds, 
while  over  five  successive  sweeps  was  less  than  1.0  x  10"3.  Intermediate  results 
were  monitored  to  insure  convergence  and  numerical  stability  at  a  specific  time 
step  during  the  course  of  the  solution.  These  intermediate  results  were  also 
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used  to  determine  convergence  toward  either  a  steady  state  or  a  steady  periodic 
solution. 

Appropriate  grid  size  and  spacing  was  chosen  to  give  acceptable  numerical 
accuracy  while  still  maintaining  reasonable  and  acceptable  computational  times. 
Generally,  the  approach  to  organizing  the  grid  layout  for  typical  dosed  cavity 
buoyant  problems  has  been  to  concentrate  on  an  increased  number  of  grids  at  the 
walls  to  accurately  describe  the  thermal  and  momentum  boundary  layers.  For  the 
two-dimensional  formulation,  a  uniform  100  x  100  grid  was  chosen  to  give  good 
coverage  without  knowing  in  advance  as  to  where  the  greatest  dependent  variable 
gradients  would  be  located  within  the  chamber  cavity.  This  grid  layout  was 
subsequ<  ntly  refined  in  regions  of  interest  to  verify  numerical  accuracy.  Refining 
the  grid  by  doubling  the  grid  number  at  the  vertical  wall,  a  high  gradient  region, 
showed  that  the  predicted  maximum  wall  boundary  layer  velodty  using  the  100  x 
100  grid  underestimated,  to  within  3  percent,  the  value  obtained  using  the  refined 
grid.  The  predicted  momentum  boundary  layer  thickness  was  less  than  5  percent 
for  boundary  layers  greater  than  2.5  cm  in  thickness  while  the  temperature  was 
less  than  1  percent  of  the  value  predicted  using  the  refined  grid. 

While  the  100  x  100  grid  in  the  two-dimensional  formulation  was  found  to  be 
satisfactory,  expansion  to  a  three-dimensional  formulation  was  found  to  be  com¬ 
putationally  intensive  and  time  consuming.  The  third-dimension  grid  number 
was  set  to  five  to  maintain  reasonable  computational  times  while  still  being  able 
to  qualitatively  describe  the  general  fiow  structure  for  the  **  plane  and  thereby 


address  some  of  the  inconsistencies  between  the  numerical  and  experimental  re¬ 
sults. 

Nusselt  number  was  not  evaluated  since  this  only  gives  information  regarding 
gross  behavior  and  is  insensitive  to  specific  information  regarding  the  flow  struc¬ 
ture  and  temperature  distribution  at  various  chamber  locations  such  as  the  core 
region  (Ostrach  and  Austin,  1988).  This  did  not  preclude  evaluation  of  the  Nus¬ 
selt  number  to  verify  numerical  convergence  and  optimum  grid  size  but  still  only 
gave  information  for  the  numerical  conditions  at  the  chamber  wall.  In  addition, 
there  is  an  error  associated  with  the  Boussinesq  approximation  (Zhong  et  al., 
1985)  in  which  there  is  as  much  as  a  20  percent  overestimation  of  the  maximum 
vertical  velocity. 

1.5  Experimental  Approach 


An  experimental  chamber  model  was  fabricated  using  1.91cm  thick  plexiglas 
with  an  interior  body  dimension  of  0.813  m  x  8.87  cm  and  0.686  m  high,  Figs. 
1.1b  and  1.2.  Baffles  of  1.5  mm  thick  plexiglas  were  inserted  across  the 
chamber  depth  and  fixed  into  the  vertical  side  walls.  The  vertical  end  walls  were 
fabricated  using  a  thin  aluminum  plate  combined  with  a  water  jacket  to  maintain 
a  constant  wall  temperature.  A  Lauda  recycling  water  bath  was  used  to  maintain 
a  constant  water  temperature  in  the  water  jacket  (±0.2°C)  Inlet  and  exhaust 
air  was  passed  through  a  series  of  laminated  porous  stainless  steel  plates  with  a 
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nominal  pore  size  of  44  microns  to  maintain  a  uniform  velocity  profile.  The  inlet 
volumetric  flow  rate  was  adjusted  and  monitored  using  a  Matheson  rotameter  and 
corrected  to  account  for  the  operational  pressure  drop.  The  chamber  pressure 
was  maintained  to  that  of  ambient  by  adjusting  the  exhaust  flow  rate. 

Field  temperature  measurements  were  made  using  iron-constantan  exposed 
junction  thermocouples  individually  calibrated  to  within  ±0.1° C.  Each  thermo¬ 
couple  was  immersed  in  a  stirred  water  bath  and  calibrated  using  a  standard 
NBS  traceable  thermometer  (Ever  Ready  Thermometer  Co.,  New  York,  N.Y.) 
(±0.05°C)  ovei  the  anticipated  range  of  operation  of  24°C  to  32°C.  An  insulated 
isothermal  thermocouple  switch  to  minimize  stray  emf  and  thermocouple  ampli¬ 
fier  were  used  in  conjunction  with  a  digital  voltmeter  for  a  sensitivity  greater 
than  the  required  ±5.1  /sV  for  ±0.1°C  measurements.  Thermocouple  voltages 
over  a  range  of  calibration  temperatures  were  obtained  and  referenced  to  an  ice 
point  consisting  of  a  large  reservoir  ice  bath  with  an  additional  reference  ther¬ 
mocouple  to  verify  shifts  in  the  reference  ice  point  temperature.  The  resulting 
calibration  was  used  to  determine  the  field  temperature  during  the  course  of  ex¬ 
perimentation.  Periodically,  the  standard  NBS  traceable  thermometer  was  also 
used  to  verify  field  temperatures  of  the  experiment.  Thermocouples  were  cho¬ 
sen  and  sized  to  minimize  flow  disturbances  and  maximize  response  time  with 
a  sheath  diameter  of  0.159  cm  and  a  bead  diameter  ranging  from  0.025  cm  to 
0.076  cm  for  a  time  constant  of  1.35  to  1.73  seconds. 

Velocity  data  was  obtained  using  a  TSI IFA-100  hot  wire  anemometer  using  a 
platinum  film  sensor  model  #1211-10.  To  effectively  use  the  hot  wire  anemome- 
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ter,  the  magnitude  of  the  buoyant-induced  velocity,  due  to  the  heating  of  the 
sensor,  in  relation  to  the  velocity  to  be  measured  was  of  interest.  Uncertainty 
errors  can  be  determined  depending  upon  the  orientation  of  the  fluid  velocity 
to  that  of  the  induced  buoyant  velocitv  as  described  by  TSI  (TSI  TB  14).  This 
is  '-vident  when  the  stream  velocity  is  opposed  to  the  censor-induced  buoyant 
velocity  and  of  the  same  order  of  magnitude.  During  this  experiment,  there  was 
a  range  of  stream  velocities  in  which  the  indicated  sensor  velocity  would  decrease 
as  the  stream  velocity  is  increased.  This  is  due  to  the  counter  acting  velocities 
summing  and  indicating  a  lower  resultant  velocity.  A  spot  calibration  of  the 
prob?  sensor  using  a  series  of  parallel  screens  at  isothermal  conditions  showed 
♦  his  velocity  to  be  on  the  order  of  2.15  x  10~J  ±  0.25  x  10~*  m/s.  Therefore 
hot  wire  anemometry  was  limited  to  the  evaluation  of  the  momentum  boundary 
layer  thickness  along  the  vertical  end  wall  with  the  edge  of  the  boundary  layer 
being  defined  as  this  baseline  velocity. 

1.0  Results  and  Discussion 


Experimental  data  from  the  experimental  model  consisted  of  temperature 
measurements  at  various  chamber  locations  at  the  eentral  plane  and  hot  wire 
anemometry  data  to  determine  the  momentum  boundary  layer  thickness.  These 
experimental  data  were  then  compared  with  results  obtained  from  the  three- 
dimensional  numerical  solutions,  also  at  the  central  plane,  to  simulate  this  same 
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experimental  model.  Two-dimensional  numerical  solutions  to  simulate  the  full- 
size  chamber  were  also  compared  to  the  same  model  experimental  data  and  the 
corresponding  three-dimensional  numerical  solutions. 

To  obtain  similitude  of  the  full-size  chamber  using  the  experimental  model, 
the  non-dimensional  parameters,  /lj/Re,  Gr/AjRe2  and  Aj/Pr-Re  were  matched 
with  the  exception  of  the  coefficient,  Aij  A\.  The  aspect  ratio  Aj  was  not  consid¬ 
ered  for  similitude  of  the  full-size  chamber  and  was  either  0.813  or  12.2  for  the 
full  size  chamber  and  experimental  chamber  respectively.  As  a  result,  the  inlet 
velocity,  Wjm,  and  temperature  difference,  6T,  between  the  inlet  and  vertical  end 
walls  remained  the  same  as  that  of  the  full-size  chamber  since  only  the  chamber 
depth  varied. 

Figures  1.3- 1.6  show  the  chamber  velocity  direction,  velocity  magnitude  and 
temperature  contours  for  the  two-dimensional  cases  #l-#4.  Figures  1.7-1.10 
show  the  temperature  contours  and  velocity  vector  and  magnitude  plots  for  the 
three-dimensional  cases  #5  and  #6.  Final  comparisons  of  non-dimensional  tem¬ 
perature  are  shown  in  Figs.  1.11  and  1.12  for  the  two-dimensional  full  size  cham¬ 
ber  case  # 4,  the  three-dimensional  cases  #5  and  #6,  and  experimental  results 
case  #7.  Momentum  boundary  layer  thickness  and  maximum  velocity  for  the 
cases  investigated  are  shown  in  Table  1.2. 

Three  characteristic  flow  structures  were  observed  for  negative,  positive  and 
zero  differences  between  the  inlet  and  vertical  wall  temperatures.  The  three- 
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dimensional  cases  exhibited  additiunai  flow  structure  characteristics  inherent  to 
the  third-dimensional  boundary  condition. 

Negative  differences  between  the  inlet  and  vertical  end  wall  temperatures 
induced  "opposed’  buoyant  forces  at  the  vertical  end  walls  as  > wn  m  Figs. 
1.3b  and  1.4  for  the  two-dimensional  case  tt  2.  The  flow  structure  for  this  case 
consists  of  a  central  downward  flowing  core  of  cooler  fluid  with  a  reflux  of  mass 
at  the  vertical  end  walls  from  the  exhaust  to  the  inlet  due  to  the  upward  vertical 
velocity  component  induced  by  the  "opposed’  buoyant  forces.  This  reflux  of 
mass  subsequently  combines  with  the  co  der  inlet  mass  at  the  inlet  plane.  The 
resulting  flow  structure  was  not  only  steady  periodic  in  nature  with  &  period  on 
the  order  of  300  seconds,  but  was  aiso  displaced  from  the  expected  symmetrical 
flow  pattern.  Placement  of  the  lower  baffles  in  relation  to  the  exhaust  plane  was 
found  to  be  the  contributing  factor  in  the  development  of  the  oscillatory  flow  and 
the  asymmetric  flow  structure.  When  these  lower  baffles  were  removed,  the  flow 
structure  became  symmetric  with  little  oscillation. 

Positive  and  incre&sinf  differences  bet—^en  the  inlet  and  end  wall  tempera¬ 
tures  resulted  in  a  increased  downward  vertical  velocity  component  at  the  vertical 
end  walls  with  mass  refluxing  from  he  exhaust  plane  and  flowing  horizontally 
along  the  upper  surface  of  the  lower  baffles.  This  is  shown  in  Figs.  1.5-1. 7  for 
the  two-dimensional  cases  9  3  ami  y  4  i.nd  the  three-dimensional  adiabutic  side 
wail  case  As  the  temperature  differences  increased,  both  a  thirning  of  the 
vertical  wall  boundary  layer  coupled  with  an  increased  vertical  velocity  compo¬ 
nent  waa  observed  as  shown  in  Table  1.2.  With  this  increase  in  the  temperature 
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differential,  the  central  vertical  mass  flow  between  the  baffles  was  diminished  as  a 
result  of  an  increased  mass  flux  to  the  vertical  end  walls  above  the  upper  baffles. 

There  are  two  distinct  regions  of  interest.  The  first  regions  of  interest  are 
the  areas  above  the  upper  baffles,  locations  where  flow  is  diverted  around  the 
baffles,  and  along  the  vertical  end  walls.  In  these  areas,  convection  begins  to 
dominate  diffusion  and  becomes  the  primary  mode  of  energy  transport.  This  is 
demonstrated  in  the  temperature  contour  plots  of  Figs.  1.5,  1.6,  and  1.7  by 
the  constant  temperature  region  and  the  distorted  contours  between  the  baffles 
themselves  and  the  chamber  walls.  The  second,  lower,  region  of  interest  is  the 
area  below  the  upper  baffles  and  above  the  lower  baffles  where  diffusion  is  the 
primary  mode  of  energy  transport  as  a  result  of  the  low  velocity  components. 
This  is  demonstrated  by  the  stratified  temperature  gradient  in  the  contour  plots 
of  Figs.  1.5,  1.6,  and  1.7.  This  second  region  varies  in  size  and  will  gradually 
diminish  as  the  inlet  and  vertical  end  wall  temperature  differences  are  positive 
and  decreasing  to  zero.  With  decreasing  positive  temperature  differences,  the 
central  vertical  mass  flow  between  the  baffles  will  increase.  This  mass  will  enter 
into  the  lower  region  and  flow  horizontally  along  the  lower  surface  of  the  upper 
baffles  toward  the  vertical  end  walls.  As  the  temperature  difference  is  further 
decreased,  approaching  zero,  the  flow  will  depict  that  of  Fig.  1.3a  in  which  the 
buoyant  forces  are  no  longer  dominant. 

Formulation  of  the  three-dimensional  case  was  considered  such  that  the  addi¬ 
tional  side  walls  for  the  third-dimension  were  either  adiabatic  case  #5  or  constant 
temperature  case  #6.  With  this  formulation,  a  solution  with  a  conduction  side 
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wall  boundary  condition,  as  was  anticipated  using  the  experimental  model,  would 
be  bordered  by  the  two  solutions  with  the  side  wall  boundary  condition  being 
either  adiabatic  or  constant  temperature.  For  constant  temperature  side  walls, 
the  temperature  was  maintained  less  than  the  inlet  temperature  and  greater  than 
the  vertical  end  wall  temperature.  Average  experimental  wall  temperatures  were 
obtained  during  the  experimental  testing  and  subsequently  used  as  the  bound¬ 
ary  condition  for  the  side  wall  temperature  in  the  three-dimensional  numerical 
formulation  of  case  #6.  For  this  case,  the  side  wall  temperature  was  maintained, 
satisfying  the  following  relation  for  comparison  to  the  non-dimensional  temper¬ 
ature,  T  : 

(Tin  ~  T,w)  =  T'(Tjn  -  Tw) 

where  T' — 1/6 

Comparison  of  the  two-dimensional  case  #4  and  the  three-dimensional  adia¬ 
batic  side  wall  case  #5  show  a  slight  variation  in  the  temperature  contour  plots 
(Figs.  1.6b  and  1.7b.)  attributed  to  the  increased  aspect  ratio,  Aj  and  the 
additional  no-slip  side  wall  boundary  condition.  Comparison  of  the  solution  of 
the  three-dimensional  constant  side  wall  temperature  case  #6  to  that  of  the 
three-dimensional  adiabatic  side  wall  case  #5  and  the  two-dimensional  case  #4 
showed  a  significant  deviation  in  both  the  flow  structure  and  temperature  profile 
as  shown  in  Figs.  1.6-1. 8.  This  variation  between  the  two  three-dimensional 
solutions  are  attributed  to  differences  in  the  side  wall  boundary  conditions. 

Additional  effects  of  the  third-dimension  can  be  further  demonstrated  from 
Figs.  1.9  and  1.10.  The  three-dimensional  adiabatic  case  #5  differs  from  the 
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two-dimensional  formulation,  case  #4,  in  that  the  no-slip  boundary  condition 
is  included  along  the  side  walls.  The  constant  temperature  side  wall  case  #6 
resulted  in  an  increased  boundary  layer  flow  that  is  characteristic  of  the  "aided” 
mixed  convection  flow  in  a  pipe  or  channel.  Case  #6  also  exhibited  flow  reversal 
in  areas  above  the  lower  baffles  and  at  the  inlet.  Both  cases  #5  and  #6  showed 
rotational  circulation  below  the  lower  baffles. 

Non-dimensional  temperature  profiles  for  two-dimensional  case  #4  and  three- 
dimensional  cases  #5  and  #6  were  compared  with  the  experimental  results  case 
#7  and  are  shown  in  Figs.  1.11  and  1.12.  In  these  figures,  non-dimensional 
center  plane  temperatures  along  the  chamber  vertical  axis  are  compared  at  lat¬ 
eral  locations  along  the  Y’-axis,  where  Y’=0  is  the  centerline  of  the  chamber 
and  Y’=l  the  location  of  the  vertical  end  wall.  The  two-dimensional  case  #4 
consistently  overestimated  the  non-dimensional  temperature  when  compared  to 
the  three-dimensional  adiabatic  case  #5.  These  differences  in  non-dimensional 
temperature  were  most  evident  for  Z/H  >  0.4  and  were  also  attributed  to  the 
increased  aspect  ratio,  Aj,  and  the  additional  no-slip  side  wall  boundary  condi¬ 
tion.  There  were  significant  differences  between  these  two  cases  (#4  and  # 5), 
case  #6  and  the  experimental  results  case  #7.  The  experimental  results,  as  ex¬ 
pected,  agree  more  with  a  conduction  side  wall  boundary  condition  as  opposed  to 
the  chosen  adiabatic  or  constant  temperature  boundary  conditions.  A  conjugate 
numerical  formulation  to  account  for  the  side  wall  boundary  condition  of  third- 
dimension  for  comparison  to  the  experimental  results  would  be  more  appropriate 
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than  the  chosen  boundary  condition,  but  would  require  a  significant  increase  in 
computational  time. 

Verification  of  the  vertical  end  wall  boundary  layer  thickness  was  performed 
using  hot  wire  anemometry  measurements.  The  numerical  error  in  boundary 
layer  thickness  was  taken  to  be  one  half  of  the  grid  size,  ±0.41  cm,  and  the 
experimented  error  ±0.32  cm.  Anemometer  measurements  were  made  by  noting 
a  change  in  the  response  from  the  zero  baseline  and  the  corresponding  location 
from  the  vertical  end  wall  determined  to  be  the  boundary  layer  thickness.  This 
boundary  layer  thickness  as  defined  was  also  determined  from  the  numerical 
solutions  of  cases  #3  -  #6  and  the  results  compared  in  Table  1.2  for  the  ver¬ 
tical  non-dimensional  locations  at  ZjH =0.44  and  ZJH— 0.67.  There  was  good 
agreement  within  error  of  both  the  experimentally  determined  boundary  layer 
thickness  and  the  numerically  predicted  boundary  layer  thickness. 

1.7  Conclusions 


Development  of  the  chamber  flow  structure  was  found  to  be  sensitive  to 
small  variations  in  the  temperature  difference  between  the  inlet  and  wall.  For 
positive  temperature  differences  greater  than  1°C,  the  entry  above  the  upper 
baffles  was  sufficient  such  as  to  provide  enough  length  for  the  flow  to  develop 
and  combine  in  the  boundary  layer  flow.  Upper  baffle  placement  in  this  instance 
appeared  to  have  minimal  effect  upon  the  development  of  the  resulting  flow 
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structure.  However,  with  gradually  decreasing  temperature  differences,  the  baffle 
placement  had  significant  impact  upon  the  flow  structure,  eventually  resulting 
in  an  asymmetric  steady  periodic  behavior.  There  were  also  distinct  regions  of 
fluid  and  energy  transport  due  either  to  diffusion  or  convection  depending  upon 
the  location  within  the  confines  of  the  chamber.  With  increasing  or  decreasing 
temperature  differentials,  these  areas  were  either  enlarged  or  diminished  as  the 
flow  structure  developed  with  temperature  being  the  driving  parameter. 

Use  of  the  laboratory  experimental  model  as  an  indicator  in  the  evaluation  of 
the  flow  structure  in  the  larger,  full-size  chamber  appeared  be  limited  due  to  the 
significant  effects  of  the  three-dimensional  geometry.  Effects  of  the  additional 
side  wall  boundary  condition  appeared  to  have  a  significant  effect  upon  the  de¬ 
velopment  of  the  flow  structure  within  the  model,  significantly  differing  from  the 
two-dimensional  numerical  solution. 

In  general,  the  observed  variations  in  development  of  the  flow  structure  from 
steady  state  to  oscillatory  behavior  are  consistent  with  laminar  buoyant  flow 
bifurcations  for  which  Prandtl  number,  geometry,  and  Rayleigh  number  are  the 
parameters  influencing  the  flow  structure  at  a  specific  bifurcation  (Yang,  1988). 
Such  appears  to  be  the  case  with  the  development  of  the  flow  structure  in  the 
chambers  investigated  here.  Variation  in  Rayleigh  number  due  to  the  difference 
between  the  inlet  and  wall  temperature  in  conjunction  with  either  the  baffle 
placement  or  chamber  aspect  ratios  provided  the  means  for  these  flow  structure 
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transitions  to  occur.  This  resulted  in  the  observed  differences  between  the  two- 
dimensional  formulation  to  simulate  the  full-size  chamber,  the  three-dimensional 
formulations  and  laboratory  experimented  model. 


Table  1.1 

PARAMETRIC  SPECIFICATIONS  FOR  NUMERICAL 
AND  EXPERIMENTAL  INVESTIGATION 
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case  # 

Approach 

T,w 

Tv 

Tin 

ST 

1 

2D  NUMERICAL 

N/A 

24°  C 

24°C 

0°C 

2 

2D  NUMERICAL 

N/A 

24°  C 

23°  C 

-1°C 

3 

2D  NUMERICAL 

N/A 

24°  C 

25°C 

1°C 

4 

2D  NUMERICAL 

N/A 

24°C 

27.3°C 

3.3°C 

5 

3D  NUMERICAL 

adiabatic 

24°  C 

27.3°C 

3.3°C 

6 

3D  NUMERICAL 

constant** 

24°  C 

27.3°C 

3.3°  C 

7 

3D  EXPERIMENTAL 

conducting 

24°C 

27.3°  C 

3.3°C 

case  # 

Ra 

Gr 

Re 

M 

Aj 

Ai 

1 

0 

0 

235 

1.19 

0.813 

1.31 

2 

(-)*6.82  x  107 

(-)9.74  x  107  235 

1.19 

0.813 

1.31 

3 

6.82  x  107 

C  74  x  107  235 

1.19 

0.813 

1.31 

4 

2.25  x  10* 

3.21  x  10*  235 

1.19 

0.813 

1.31 

5 

5.84  x  105 

8.35  x  10*  32.3 

1.19 

12.2 

0.180 

6 

5.84  x  10s 

8.35  x  10s  32.3 

1.19 

12.2 

0.180 

7 

5.84  x  10s 

8.35  x  105  32.3 

1.19 

12.2 

0.180 

*  negative  sign  signifies  "opposed”  and  positive  sign  signifies  "aided”  buoyant 
force  at  vertical  "end  walls” 


’*  T,w=(TlN-6T/6)ot  T- 1/6 
Pr=0.7 


W[ ft  =  4.1  x  KT1  m/s 
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Table  1.2 

BOUNDARY  LAYER  THICKNESS  AND  MAXIMUM  VELOCITY 


case  # 

Vnax 

6 *,  cm 

characteristic 

Z 1 H =0A4 

Z/H= 0.67 

error,  cm 

1 

2.81  x  10"J 

N/A 

N/A 

N/A 

steady  state 

2 

7.03  x  10"3 

N/A 

N/A 

N/A 

steady  periodic 

3 

6.89  x  10"3 

2.38 

2.84 

±0.41 

steady  state 

4 

1.03  x  10" 1 

1.98 

2.32 

±0.41 

steady  state 

5 

1.02  x  10" 1 

1.93 

2.31 

±0.41 

steady  state 

6 

Sj 

1.16  x  10'1 

2.09 

2.58 

±0.41 

steady  state 

7 

not  measured 

1.91 

2.54 

±0.32 

steady  state 

•6  determined  to  be  the  location  for  a  velocity  equivalent  to  the  minimum 
response  of  the  anemometer  sensor  («  2.15  x  10"3  m/s) 
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Fig.  1.1  Chambr-  if  ■Miietry  ipecification*  (a)  two-dimensional  geometry 
(b)  three-diraen»iond  laboratory  model  geometry 


Fig.  1.5  Two-dimensional  results  for  case  #3  (a)  chamber  velocity  direc¬ 
tion  and  magnitude  profile  (b)  temperature  contour 
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Fig.  1.6  Two  dimensional  result*  for  ca«r  ff  \  (a)  chamber  velocity  direc¬ 


tion  and  magnitude  profile  (b)  temperature  contour 
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Fig.  1.8  Three-dimensional  results  for  case  #6  at  center  plane  (a)  chamber 


velocity  direction  and  magnitude  profile  (b)  temperature  contour 


rig.  1.10  Velocity  vector  plot  acrosa  depth  of  three-dimenrional  numerical 
t ulation  at  lateral  distance  .251  aa  viewed  from  the  origin  (a)  caae  #5 
caae  #6 


Chapter  2 


MIXED  CONVECTION  ANALYSIS  IN  LARGE  BAFFLED 
RECTANGULAR  CHAMBERS  WITH  INTERNAL  HEAT  SOURCES 


2.1  Summary 


An  investigation  to  determine  the  effects  of  mixed  convection  on  the  flow 
structure  of  large  baffled  chambers  with  internal  heat  sources  has  been  com¬ 
pleted.  Baffles  and  internal  heat  sources  were  placed  symmetrically  about  the 
vertical  axis  while  maintaining  a  constant  and  uniform  inlet  velocity,  inlet  tem¬ 
perature,  and  wall  temperature  of  4.1  x  10“*  m/s,  27*C  and  24*C  respectively. 
Time  dependent  experimental  and  two-dimensional  numerical  results  were  com¬ 
pared  at  various  chamber  locations  to  demonstrate  the  transition  of  the  flow 
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structure  through  a  sequence  of  bifurcations  with  increasing  internal  heat  source 
temperature  up  to  32°C. 

Phase-space  and  frequency  plots  cf  temperature  and  velocity  for  the  numer¬ 
ical  solutions  showed  that  there  were  characteristic  harmonics,  amplitude,  and 
phase-space  shapes  at  the  various  chamber  locations.  The  three-dimensional  ex¬ 
perimental  laboratory  model  demonstrated  a  transition  of  the  flow  structure  from 
one  that  is  steady  state  to  periodic,  aperiodic  and  Anally  chaotic  with  increasing 
internal  heat  source  temperature.  Three-dimensional  experimental  results  also 
exhibited  unique  harmonics  and  amplitude  of  temperature  and  velocity  at  vari¬ 
ous  chamber  locations.  Although  differing  from  the  two-dimensional  numerical 
results  at  the  same  chamber  locations,  the  harmonics  observed  from  the  three- 
dimensional  experimental  model  showed  similar  frequency  behavior  and  trends 
to  that  of  the  two-dime;  onal  numerical  restilts. 

2.2  Introduction 


Large  baffled  rectangular  chambers  with  low  inlet  velocities  have  traditionally 
been  used  to  study  a  variety  of  chemical  atmospheres  generated  for  health  effects 
and  biological  studies,  chemical  reaction  and  chemical  species  formation  stud¬ 
ies,  and  aerosol  characterization  studies.  Generally,  they  have  a  cross-sectional 
diameter  greater  than  .5  meter  with  an  inlet  velocity  set  to  maintain  10-15  cham¬ 
ber  volume  changes  per  hour.  The  inlet  and  exhaust  may  be  mounted  vertically 
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or  hoiizontaHy  with  an  array  of  internal  baffles  included  for  atmosphere  disper¬ 
sion  purposes.  Additional  capabilities  include  cage  units  for  housing  laboratory 
animals  for  a  variety  of  biological  research  applications. 

There  have  been  many  investigations  concerning  the  Rayleigh- Benard  flow 
in  confined  and  partially  open  enclosures  with  and  without  baffles  or  dividers 
and  internal  or  wall  heat  sources  (Bajorek  and  Lloyd,  1982;  Markatos  and  Ma- 
lin,  1982;  Powe  and  Warrington,  1983;  Nansteel  and  Greif,  1981;  Sparrow,  et 
al.  ,  1984;  Warrington  and  Powe,  1985;  November  and  Nansteel,  1987;  Lee  and 
Goldstein,  1988;  Ostrach  and  Austin,  1988).  Such  investigations  have  made  use 
of  both  experimental  and  numerical  results  for  information  with  regard  to  the 
temperature  profile  and  flow  structure  in  various  closed  cavity  geometries  with 
differing  magnitudes  of  Grashof,  Prandtl,  and  Rayleigh  numbers. 

Yang  (1988)  discussed  the  buoyant  enclosure  problem  as  a  dissipative  dy¬ 
namic  system  whose  flow  structure  development  is  one  of  a  series  of  bifurcations 
brought  about  by  the  driving  parameters  of  Rayleigh  number,  Prandtl  num¬ 
ber,  and  enclosure  geometry.  Development  of  a  flow  structure  from  one  that 
is  steady  state,  periodic  in  time,  aperiodic,  and  finally  chaotic  are  examples  of 
specific  bifurcations  that  may  be  expected  with  variation  in  Rayleigh  number, 
Prandtl  number,  and  altered  geometry.  Yang  (1988)  further  discussed  the  use 
of  an  attractor,  a  concept  used  in  nonlinear  dynamics,  to  describe  the  temporal 
asymptotic  behavior  of  a  trajectory  in  the  phase  space  corresponding  to  a  specific 
bifurcation. 
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Initial  investigations  of  these  chambers  incorporated  flow  visualization  tech¬ 
niques  using  smoke  and  dyes  to  optimize  the  chamber  geometry  so  that  well- 
dispersed  atmospheres  could  be  obtained  (Carpenter  and  Beethe,  1981;  Moss, 
1981).  Quantitative  approaches  to  determine  the  mean  residence  time  of  a  cham¬ 
ber,  the  rate  of  dispersion  within  a  chamber  of  a  tracer  gas,  and  the  distribution 
within  a  chamber  using  vapors,  droplets,  and  solid  particles  have  also  been  done 
(Hemenway  et  al.,  1982;  Beethe  et  al.,  1979;  Moss,  1982;  Yell  et  al.,  1986).  Most 
recently,  Yerkes  et  al.  (1989)  investigated  the  mixed  convection  in  large  baf¬ 
fled  rectangular  chambers  and  the  characteristic  flow  structure  development  for 
"aided”  and  "opposed”  buoyant  forces.  Development  of  the  chamber  flow  struc¬ 
ture  was  found  to  be  predominantly  buoyant  in  nature  and  sensitive  to  smell 
variations  in  the  temperature  difference  between  the  inlet  and  wall  varying  from 
steady  state  to  asymmetric  oscillatory  behavior.  These  investigations,  however, 
have  not  addressed  the  effects  of  mixed  convection  with  internal  biological  heat 
sources  on  the  transport  phenomenon  of  heat  and  mass  transfer  within  ther 
chambers.  Information  into  the  diffusion  and  convection  of  mass,  enthalpy,  and 
momentum  within  these  chambers  with  the  biological  heat  sources  are  of  inter¬ 
est  to  better  understand  their  operation  and  provide  a  basis  for  improvement, 
modification  and  future  development. 

A  typical  chamber  that  may  be  used  for  the  purposes  described  consists 
of  a  conical  inlet  intended  to  uniformly  distribute  the  atmosphere  of  interest 
aeross  the  inlet  plane.  The  exhaust  consists  of  a  manifold  to  achieve  uniform 
flow  chsuacteristics  across  the  exhaust  plane.  Inlet  and  exhaust  flow  rates  are 
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controlled  independently  to  maintain  a  slight  negative  chamber  pressure  with 
respect  to  the  ambient  pressure  (2.5  -  5.1  cm  H2O).  Biological  heat  sources  may 
consist  of  a  variety  of  laboratory  animal  species  and  are  housed  in  uniformly 
spaced  cage  units  located  directly  over  the  baffles. 

This  investigation  considers  a  simplified  version  of  the  aforementioned  cham¬ 
ber.  Of  interest  is  the  transient  behavior  and  development  of  the  flow  structure 
due  to  a  temperature  gradient  between  the  biological  heat  sources,  inlet  tem¬ 
perature  and  wall  temperature.  Results  from  a  time  dependent  two-dimensional 
numerical  formulation  were  compared  to  experimental  results  from  a  scaled-down 
laboratory  experimental  model  intended  to  be  used  as  an  indicator  for  use  in  the 
analysis  of  the  larger,  full-size  chamber.  The  two-dimensional  numerical  formu¬ 
lation  corresponded  to  the  larger,  full-size  chamber  differing  from  that  of  the 
experimental  laboratory  model  in  that  the  depth  is  1.0  m  compared  to  6.67  cm 
for  the  model.  The  vertical  end  wall  temperature,  T»,  was  assumed  to  be  con¬ 
stant  and  equal  to  24°C  due  to  the  thin  wall  construction.  The  inlet  temperature 
was  alro  assumed  to  be  constant  across  the  inlet  plane  and  equal  to  27°C. 

Biological  heat  sources  were  treated  as  a  constant  temperature  surface  at 
32°C  and  of  a  size  and  mass  that  would  be  appropriate  for  a  200  gm  to  300  gm 
rat  in  a  confined  enclosure  (Mauderly,  1986).  It  was  assumed  that  physiolog¬ 
ical  control  of  the  body  temperature  is  sufficient  to  maintain  a  constant  body 
temperature  in  the  low  air  velocities  within  the  confines  of  the  chamber. 
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Typical  values  of  Grashof,  Rayleigh  and  Reynolds  numbers  for  this  study 
were  Gr  =  2.02  x  10s,  Ra  =  1.42  x  10s  and  Re  =  32  and  235  using  the  inlet 
hydraulic  diameter,  D as  the  characteristic  length  (the  two-dimensional  formu¬ 
lation  assumes  chamber  depth  of  one  meter  while  the  experimental  model  has 
a  depth  of  6.67  cm  thereby  resulting  in  different  Reynolds  numbers  due  to  the 
differing  hydraulic  diameters).  The  difference  in  temperature  was  taken  to  be  the 
maximum  difference  between  the  vertical  end  wall  temperature  and  the  biological 
heat  source  temperature,  ST  =  8.0°C,  with  an  inlet  velocity  of  4.1  x  10~3  m/s. 
With  these  values  for  Grashof,  Rayleigh  and  Reynolds  numbers,  the  development 
of  the  flow  structure  is  predominately  buoyant  in  nature  with  the  buoyant  forces 
dominating  the  inertia  forces. 

2.3  Analysis 


Two-dimensional  numerical  solutions  of  the  incompressible  time-dependent 
laminar  Navier-Stokes  equations  were  calculated  in  accordance  with  the  two- 
dimensional  problem  shown  in  Fig.  2.1a.  For  this  investigation,  only  the 
dimensional  form  of  the  governing  equations  were  solved,  thereby  circumventing 
the  immediate  need  for  scaling  considerations  of  the  non-dimensional  parameters 
as  is  typically  required  for  buoyant  flow  conditions  (Ostrach  and  Austin,  1988). 
The  Boussinesq  approximation  was  used  in  the  vertical  direction  to  account  for 
buoyancy  effects.  Using  this  approximation  the  vertical  body  force  becomes: 
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Bt  =  pg[l-HT-TlN)} 


where 


for  an  ideal  gas. 


Fcr  this  investigation,  the  reference  Talues  for  all  of  the  gas  properties  were 
taken  to  be  those  at  the  inlet  plane. 

The  chamber  inlet  cone  was  assumed  to  provide  a  uniform  velocity  distri¬ 
bution  across  the  inlet  plane  with  the  exhaust  manifold  providing  a  velocity 
distribution  resulting  from  a  constant  exhaust  pressure  across  the  outlet  plane. 
Viscous  dissipation  and  pressure  work  were  assumed  to  be  negligible  due  to  the 
low  buoyancy-induced  velocities.  With  these  assumptions,  the  two-dimensional 
governing  equations  become: 


Conservation  of  mass: 


dv  dw 


(2.1) 


Conservation  of  momentum: 


(  dv  dv  8v\  6p*  ( d3v  d3v\ 

'  U  +  v  a;  +  "  a;  J=  -  97  +  *  (  v  +  a?  J 


(2.2) 


p 


dw  dw  dw\  dp+  ( d2w  d2w\ 

Ji  +  va;  +  wTz)  =  "  aT  +  'W+a?j 


-rpg3{T  -  T/.v) 


(2.3) 


Conservation  of  energy: 

(dT  dT  dT\  k  (  d2T  d2T\ 

'(dT  -  c~[a?+e?) 


(2.4) 


The  boundary  conditions  were  specified  with  the  assumption  that  there  is  a 
uniform  inlet  velocity  profile  and  the  vertical  end  wall  temperature  is  constant. 
Only  the  magnitude  of  the  inlet  velocity  is  specified,  lhe  magnitude  of  the  outlet 
velocity  is  such  as  to  satisfy  the  conservation  of  mass  with  a  constant  reference 
pressure  across  the  outlet  plane. 

The  outlet  boundary  condition  for  temperature  was  considered  to  be  locally 
parabolic  such  that  the  Pedet  number  is  sufficiently  large  so  as  to  exhibit  lo¬ 
cal  one-way  behavior  in  the  axial  direction.  The  baffles  were  treated  as  being 
infinitesimally  thick  with  a  no-slip  boundary  condition.  The  baffle  temperature 
was  specified  to  be  that  of  the  local  fluid  temperature. 


Boundary  conditions: 


Inlet:  t>0,  v  ~  0,  w  =  Win  ~  4.1  x  10''*  m/s,  T  =  Tin  =  27°C 


Vertical  end  wall:  t>0,  v  =  0,  w  -  0,  T  =  T„  =  24°C 
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Baffle:  t>0,  v  =  0,  w  =  0 


Heat  sources:  t=>0,  T  =  32°C 


Outlet:  t>0,  v  —  0,  p+  =  0 


For  this  investigation,  two  restart  fields  were  chosen  for  the  initial  conditions 
of  the  field  dependent  variables.  The  method  for  obtaining  the  restart  fields  are 
discussed  in  the  following  section. 

2.4  Numerical  Scheme 


The  numerical  scheme  used  is  a  finite-difference  iterative  method  of  solution 
using  a  control- volume  approach  as  developed  by  Spalding  et  al.  (1980).  The 
two-dimensional  numerical  model  was  formulated  using  a  fully  implicit  finite- 
difference  scheme  of  the  time  dependent,  laminar  Navier-Stokes  equations.  The 
"SIMPLEST”  (Spalding  et  al.,  1980)  method  of  solution  for  the  momentum 
equations  was  used  with  the  hybrid  differencing  formulation.  When  a  cell  Peclet 
number  is  within  the  range  -2  to  2,  a  central-difference  scheme  is  used  and  when 
the  cell  Peclet  number  is  outside  this  range  the  upwind  differencing  scheme  is 
used. 
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The  outlet  boundary  condition  for  temperature  was  considered  to  be  locally 
parabolic  such  that  the  Peclet  number  is  sufficiently  large  so  as  to  exhibit  local 
one-way  behavior  in  the  axial  direction  (Patankar,  1980).  This  approach  to 
specifying  the  outlet  boundary  condition  for  a  dependent  variable  is  generally 
accepted  as  the  most  available  approach  for  internal  flow  problems  in  which 
the  fluid  leaves  the  calculation  domain  and  the  dependent  variable  of  interest  is 
unknown. 

The  baffles  were  treated  as  being  infinitesimally  thick  with  a  no-slip  bound¬ 
ary  condition  by  setting  the  cell  wall  velocity  to  zero  over  the  region  of  the  baffle 
location.  The  baffle  temperature  was  specified  to  be  that  of  the  local  fluid  tem¬ 
perature.  Heat  sources  were  ’’blocked”  by  setting  the  cell  velocity  components 
within  the  heat  source  region  to  zero  and  fixing  to  a  constant  temperature. 

The  solution  sequence  involved  solving  for  the  velocity  and  temperature  fields 
from  an  assumed  pressure  field  at  a  specific  time  step.  The  pressure  field  was 
then  subsequently  updated  using  these  velocity  and  temperature  fields  in  the 
pressure-correction  equation  such  that  the  conservation  of  mass  was  satisfied. 
This  iterative  sweep  process  was  then  repeated  using  a  slab-by-slab  method  at 
each  time  step  until  convergence  criteria  were  satisfied.  Convergence  character¬ 
istics  at  each  time  step  were  maintained  by  reducing  the  time  step  where  appro¬ 
priate.  Conditions  for  convergence  to  a  solution,  either  at  a  specific  time  step  or 
to  a  steady  state  solution,  was  based  upon  a  dependent  variable,  <f>,  varying  less 
than  a  predetermined  change  of  the  magnitude  fraction,  e^,  between  successive 
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sweeps  or  time  steps  where: 


I  <t>s  -  4i+i  I  I4><e+ 

The  values  of  were  varied  and  it  was  required  that  the  magnitude  of  the 
dependent  variable  be  greater  than  the  roundoff  error  or  truncation  error  such 
as  would  happen  when  velocities  approach  zero.  Typically,  was  less  than 
1.0  x  10~5  for  five  successive  sweeps.  Intermediate  results  were  monitored  to 
insure  convergence  and  numerical  stability  at  a  specific  time  step  during  the 
course  of  the  solution.  These  intermediate  results  were  also  used  to  determine 
convergence  toward  either  a  steady  state  or  a  steady  periodic  solution. 

Appropriate  grid  size  and  spacing  and  numerical  time  step  was  chosen  to  give 
acceptable  numerical  accuracy  while  still  maintaining  reasonable  and  acceptable 
computational  times.  Due  to  the  complexity  of  the  flow  structure,  the  errors  as  a 
result  of  grid  size  and  time  step  were  determined  first,  by  comparing  the  shape  of 
the  flow  structure  for  similarity  with  a  previous  solution  at  a  larger  time  step  or 
grid  size.  Secondly,  random  locations  throughout  the  chamber  were  selected  and 
the  average  error  compared  between  the  similar  solutions  such  that  halving  the 
time  step  resulted  in  a  less  than  a  three  percent  error  in  the  dependent  variable 
over  the  time  span  of  interest  and  decreasing  the  grid  size  by  ten  percent  resulted 
in  less  than  a  one  percent  error  in  the  dependent  variable.  For  this  problem,  a 
uniform  100  x  100  grid  was  chosen  to  give  good  coverage  with  a  time  step  of  0.010 
seconds. 


51 


Initially  the  solution  was  started  from  a  zero  velocity  and  constant  tem¬ 
perature  field.  The  time  step  was  gradually  reduced  once  the  solution  became 
oscillatory  in  nature  either  due  to  numerical  perturbations  or  naturally  occur¬ 
ring  oscillations  due  to  the  physics  of  the  problem.  Two  restart  fields  for  velocity, 
temperature  and  pressure  were  then  chosen  from  these  solutions  and  the  time  de¬ 
pendence  behavior  observed  over  a  span  of  forty  seconds  using  the  final  reduced 
time  step  of  0.010  seconds. 

With  the  development  of  the  relatively  complex  transient  flow  structure,  both 
decreasing  the  grid  size  and  time  step  would  reveal  both  additional  wavelengths 
and  frequencies  that  would  otherwise  not  be  observed  with  a  larger  grid  size  and 
time  step.  In  addition,  error  propagation  in  the  solution  with  the  progression 
of  time  would  tend  to  numerically  oscillate  the  solution  with  some  numerically 
induced  low  frequency  harmonics.  Increasing  the  grid  number  and  reducing  time 
step  is  the  appropriate  response  to  correcting  and  revealing  these  additional  wave¬ 
lengths  and  harmonics  and  decreasing  time  step  related  propagation  errors.  Due 
to  the  physical  size  of  the  problem,  however,  this  would  significantly  increase  the 
computational  time  in  solving  the  problem.  Realizing  these  numerical  limita¬ 
tions  in  obtaining  solutions  to  a  random  and  marginally  stable  problem  such  as 
presented  here  requires  such  numerical  solutions  to  be  used  with  care. 
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2.5  Experimental  Approach 


An  experimental  chamber  model  was  fabricated  using  1.91cm  thick  plexi- 
glas  with  an  interior  body  dimension  of  0.813  m  x  6.67  cm  and  0.686  m  high, 
Figs.  2.1b  and  2.2.  Baffles  of  1.5  mm  thick  plexiglas  were  inserted  across  the 
chamber  depth  and  fixed  into  the  vertical  side  walls.  The  vertical  end  walls  were 
fabricated  using  a  thin  aluminum  plate  combined  with  a  water  jacket  to  maintain 
a  constant  wall  temperature.  A  Lauda  recycling  water  bath  was  used  to  main¬ 
tain  a  constant  water  temperature  in  the  water  jacket  (±0.2°C).  Heat  sources 
were  simulated  using  flask  heating  mantles  (Glas-col)  located  symmetrically  ap¬ 
proximately  1.4  cm  above  the  baffles  and  maintained  at  constant  temperature 
(±0.5°C).  Inlet  and  exhaust  air  was  passed  through  a  series  of  laminated  porous 
stainless  steel  plates  with  a  nominal  pore  size  of  44  microns  to  maintain  a  uniform 
velocity  profile.  The  inlet  volumetric  flow  rate  was  adjusted  and  monitored  us¬ 
ing  a  Matheson  rotameter  and  corrected  to  account  for  the  operational  pressure 
drop.  The  chamber  pressure  was  maintained  to  that  of  ambient  by  adjusting  the 
exhaust  flow  rate. 

Field  temperature  measurements  were  made  using  iron-constantan  exposed 
junction  thermocouples  individually  calibrated  to  within  ±  0.1°C.  Each  thermo¬ 
couple  was  immersed  in  a  stirred  water  bath  and  calibrated  using  a  standard 
NBS  traceable  thermometer  (Ever  Ready  Thermometer  Co.,  New  York,  N.Y.) 
(0.05°C)  over  the  anticipated  range  of  operation  of  24°C  to  32°C.  An  insulated 
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isothermal  thermocouple  switch  to  minimize  stray  emf  and  thermocouple  ampli¬ 
fier  were  used  in  conjunction  with  a  digital  voltmeter  for  a  sensitivity  greater 
than  the  required  ±5.1  /iV  for  ±0.1°C  measurements.  Thermocouple  voltages 
over  a  range  of  calibration  temperatures  were  obtained  and  referenced  to  an  ice 
point  consisting  of  a  large  reservoir  ice  bath  with  an  additional  reference  ther¬ 
mocouple  to  verify  shifts  in  the  reference  ice  point  temperature.  The  resulting 
calibration  was  used  to  determine  the  field  temperature  during  the  course  of  ex¬ 
perimentation.  Periodically,  the  standard  NBS  traceable  thermometer  was  also 
used  to  verify  field  temperatures  of  the  experiment. 

Thermocouples  were  chosen  and  sized  to  minimize  flow  disturbances  and 
maximize  response  time  with  a  sheath  diameter  of  0.159  cm  and  a  bead  diameter 
ranging  from  0.025  cm  to  0.070  cm  for  a  time  constant  of  1.35  to  1.73  seconds. 
Assuming  a  maximum  time  constant  of  1.7  seconds,  the  frequency  response  for 
measuring  a  .017  Hz  (1  cycle/60  s)  signal  would  result  in  a  .985  attenuation  with 
a  -10°  phase  shift. 

Velocity  data  was  obtained  using  a  TSI  IFA-100  hot  wire  anemometer  using 
a  platinum  film  sensor  model  #1211-10  with  a  frequency  response  of  300  kHz.  To 
effectively  use  the  hot  wire  anemometer,  the  magnitude  of  the  buoyant-induced 
velocity,  due  to  the  heating  of  the  sensor,  in  relation  to  the  velocity  to  be  mea¬ 
sured  was  of  interest.  Uncertainty  errors  can  be  determined  depending  upon 
the  orientation  of  the  fluid  velocity  to  that  of  the  induced  buoyant  velocity  aa 
described  by  TSI  (TSI  TB  14).  This  is  evident  when  the  stream  velocity  is  either 
opposed  to  or  in  the  same  direction  and  of  the  same  order  of  magnitude  as  the 
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sensor-induced  buoyant  velocity.  During  this  experiment,  the  range  of  stream  ve¬ 
locities  encompassed  the  sensor-induced  velocity,  limiting  the  use  of  the  hot  wire 
anemometer  from  quantitatively  evaluating  the  velocity  magnitude  to  evaluating 
the  frequency  behavior  of  the  velocity  at  various  chamber  locations. 

The  laboratory  experimental  chamber  was  operated  in  such  a  manner  so  as 
to  observe  the  development  of  the  transient  behavior  of  both  temperature  and 
velocity  at  six  chamber  locations  (Fig.  2.1).  Initially  the  chamber  was  started 
from  a  steady  state  initial  condition  without  heat  sources  and  with  an  inlet  ve¬ 
locity,  Wjff,  of  4.1  x  10" 1  m/s.  The  heat  sources  were  gradually  brought  to  32°C 
while  maintaining  the  inlet  temperature,  T/jy ,  and  vertical  wall  temperature,  Tw, 
at  27°  C  and  24°  C  respectively.  Temperature  and  velocity  data  were  obtained  at 
the  six  chamber  locations  until  the  chamber  appeared  to  maintain  a  steady  pe¬ 
riodic  state  in  which  the  magnitude  and  oscillatory  behavior  of  velocity  and 
temperature  remained  constant.  It  was  found  that  simultaneous  measurements 
of  velocity  and  temperature  at  the  same  chamber  location  induced  temperature 
measurement  errors  due  to  the  hot  wire  anemometer  probe,  therefore,  measure¬ 
ment  .  of  temperature  and  velocity  were  obtained  at  differing  chamber  locations. 
Further  temperature  sampling  artifacts  due  to  the  hot  wire  anemometer  probe 
were  verified  to  be  negligible  by  relocating  the  hot  wire  anemometer  probe  to 
another  chamber  location  while  fixing  the  thermocouple  chamber  location. 


2.6  Results  and  Discussion 
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Experimental  data  from  the  laboratory  model  consisted  of  time-dependent 
temperature  and  hot  wire  anemometry  measurements  at  various  chamber  loca¬ 
tions  about  the  central  plane  to  determine  the  transient  behavior  of  temperature 
and  velocity.  These  experimental  data  were  then  compared  with  results  obtained 
from  the  transient  two-dimensioned  numerical  solutions. 

Figures  2.3  and  2.4  show  typical  two-dimensional  numerical  solutions  for 
velocity  and  temperature  at  a  time  t  =  t0  and  t=to+30$  respectively.  Figure  2.5 
shows  the  two-dimensional  numerical  time-dependent  variation  of  temperature 
and  velocity  for  six  chamber  locations,  ((.501,,  .31 H),  (.50 L,  .82 H),  (.70 L,  .31  //), 
(.70 L,  .62 H),  (.90/,,  .11//),  (.90 L,  .72//)].  These  numerical  solutions  were  initi¬ 
ated,  as  discussed  in  the  numerical  scheme,  from  the  two  arbitrarily  chosen  initial 
restart  fields  spanning  a  forty  second  time  frame.  Shown  in  Fig.  2.6  are  the 
numerical  solution  phase-space  shapes  which  are  plotted  as  the  rate  of  change  of 
the  dependent  variable,  velocity  and  temperature,  versus  the  dependent  variable 
over  time.  Experimentally  obtained  temperature  and  velocity  information  with 
increasing  heat  source  temperature  at  chamber  location  (.70 L,  .31 H)  are  shown 
in  Fig.  2.7  and  at  the  six  chamber  locations,  Fig.  2.8.  In  these  figures,  48 
second  snapshots  were  taken  at  10  minute  intervals  for  120  minutes  to  demon¬ 
strate  the  transient  behavior  of  the  chamber  as  the  heat  source  temperature  is 
increased.  These  experimental  data  are  also  compared  to  the  time-dependent 
velocity  obtained  from  the  numerical  results  which  are  also  shown  in  Fig.  2.8. 

Development  of  the  flow  structure,  Figs.  2.3  and  2.4,  was  found  to  be  time- 
dependent  in  nature  with  the  formation  of  four  Benard  cells  in  the  upper  portion 
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of  the  chamber  above  the  upper  baffles  and  heat  sources.  Over  the  course  of 
time,  these  ceils  dissipated  and  reformed  resulting  in  an  oscillating  core  flow 
between  the  two  upper  baffles.  A  high  velocity  wall  boundary  layer  flow  was 
maintained  on  the  vertical  end  wall  increasing  in  thickness  as  the  core  flow  dis¬ 
sipated.  Benard  cells  also  formed  between  the  upper  and  lower  sets  of  baffles 
although  the  formation  of  these  cells  appeared  to  be  influenced  by  the  core  flow. 
There  was  a  reflux  of  mass  between  the  vertical  walls  and  upper  baffles  that  was 
also  influenced  by  the  oscillating  core  flow  which  in  turn  influenced  the  forma¬ 
tion  of  the  Benard  cells  above  the  upper  baffles  and  heat  sources.  This  behavior 
resulted  in  a  marginally  stable  and  oscillatory  flow  structure. 

Two-dimensional  numerical  results  showed  the  temperature  and  velocity  to 
behave  in  an  oscillatory  nature  with  some  incommensurable  harmonics  at  the 
six  chamber  locations  investigated,  Fig.  2.5.  The  magnitude  and  frequency  of 
the  oscillations  were  shown  to  be  dependent  upon  the  chamber  location  and  the 
proximity  of  that  chamber  location  to  a  Benard  cell.  The  greatest  oscillatory 
behavior  occurred  at  the  chamber  centerline  between  the  upper  baffles,  with  the 
core  flow  developing  and  subsequently  dissipating  as  a  result  of  the  Benard  cell 
formation  above  the  upper  baffles  and  heat  sources.  Temperature  and  velocity 
tended  to  be  inversely  proportional  as  would  be  expected  in  flow  structures  driven 
by  buoyant  forces  with  a  fluid  Prandtl  number  approaching  one  for  which  the 
rate  that  the  velocity  profile  develops  is  equivalent  to  the  rate  that  the  thermal 
profile  develops,  Figs.  2.5a,  2.5b,  and  2.5d.  However,  Figs.  2.5c,  2.5e,  and  2.5f 
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did  not  show  this  behavior  possibly  due  to  the  proximity  of  the  chamber  location 
to  the  movement  of  the  low  velocity  core  of  a  Benard  cell. 

Phase-space  or  attractor  plots  of  temperature  and  velocity,  shown  in  Fig. 
2.6,  exhibited  unique  shapes  at  the  six  chamber  locations  resulting  from  a  com¬ 
bination  of  steady  state  behavior  which  is  represented  by  a  fixed  point  in  the 
phase-space,  periodic  behavior  which  is  represented  by  a  limit  cycle  in  the  phase- 
space  and  aperiodic  behavior  which  is  represented  by  a  torus  in  the  phase-space 
of  temperature  and  velocity. 

The  experimental  results  demonstrated  the  transition  of  the  flow  structure 
through  specific  bifurcations,  with  the  heat  source  temperature  being  the  driving 
parameter  for  these  transitions  to  occur.  Figure  2.7  illustrates  this  behavior  at 
chamber  location  (.70 L,  .31//).  With  an  increasing  heat  source  temperature, 
the  velocity  was  shown  to  pass  through  transitions  from  steady  state,  periodic, 
aperiodic,  and  finally  chaotic  or  random  oscillatory  behavior  with  high  frequency 
incommensurable  harmonics.  Figure  2.8  compares  the  experimentally  obtained 
temperature  and  velocity  harmonics  et  the  ».‘x  chamber  locations  to  the  velocity 
harmonics  obtained  from  the  two-dimensional  numerical  results.  The  oscillatory 
behavior  of  temperature  was  significantly  attenuated  due  to  the  frequency  re¬ 
sponse  of  the  thermocouples  and  therefore  difficult  to  resolve  and  compare  with 
velocity.  However,  the  velocity  behavior  demonstrated  unique  differences  in  har¬ 
monics  between  the  the  six  chamber  locations. 
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Experimented  results  showed  chamber  location  (.70 Z,  .62 #),  Fig.  2.8e,  to 
have  the  greatest  variation  in  velocity  while  chamber  location  (.50Z,  .82 #),  Fig. 
2.8d,  the  least.  This  is  opposite  to  the  numerical  results  which  showed  the  great¬ 
est  variation  in  velocity  to  be  at  chamber  locations  (.50 Z,  .31#)  and  (.50Z,  .82 #) 
and  the  least  variation  in  velocity  at  chamber  location  (.70Z,  .62 #).  Experimen¬ 
tal  results,  however,  did  show  that  chamber  location  (.50Z,  .31#),  Fig.  2.8a,  to 
have  the  greatest  variation  in  velocity  second  to  chamber  location  (.70Z,  .62#), 
Fig.  2.8e,  agreeing  with  the  numerical  results  as  to  being  one  of  the  two  chamber 
locations  with  the  greatest  velocity  variation. 

The  experimental  results  showed  unique  transient  behavior  at  each  of  the 
six  chamber  locations.  Transition  of  velocity  through  a  series  of  bifurca  'ons 
showed  a  unique  behavior  in  both  amplitude  and  harmonics  up  to  and  including 
the  steady  oscillatory  behavior.  Experimental  results  also  indicated  a  marginally 
stable  flow  structure  that  appeared  to  be  easily  perturbed  by  heat  source  tem¬ 
perature.  Comparison  of  the  numerically  obtained  velocity  results  to  the  exper¬ 
imentally  obtained  veiocity  data  show  good  agreement  in  harmonic  behavior  for 
chamber  locations  (.50Z,  .31#),  (.90Z,  .11#)  and  (.90Z,  .72#)  in  Figs.  2.5a, 
2.5c  and  2.5f,  respectively.  The  numerical  results  appeared  to  have  less  high  fre¬ 
quency  harmonics  than  the  experimental  results  due  in  part  to  possible  numerical 
damping.  Comparison  of  the  numerical  and  experimental  results  for  chamber  lo¬ 
cations  (.50Z,  .82#),  (.70 Z,  .31 II)  and  (.70Z,  .62#)  in  Figs.  2.5d,  2.5b  and  2.5e, 
respectively,  showed  poor  agreement  in  the  harmonic  behavior  of  velocity. 
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Verification  of  the  high  velocity  boundary  layer  region  was  performed  by 
traversing  the  hot  wire  anemometer  probe  toward  the  vertical  end  wall.  The 
location  of  this  high  velocity  region  was  determined  when  there  was  an  observable 
increase  in  the  indicated  velocity  from  that  of  a  baseline  velocity  outside  the 
confines  of  the  high  velocity  region.  Figure  2.9  shows  the  results  of  these 
measurements  confirming  the  presence  of  a  high  velocity  region  between  the 
vertical  end  wall  and  0.64±0.16  to  1.27n0.16  cm  with  a  vertical  location  of  0.67 H . 

Efforts  to  utilize  numerical  solutions  for  a  time-dependent  comparison  of  ex¬ 
perimental  results  in  a  random  and  marginally  stable  flow  structure  resulted 
in  the  comparison  of  numerical  and  experimental  techniques  that  are  each  lim¬ 
ited.  The  numerical  grid  size  and  time  step  increment  appropriate  to  capture 
wave  lengths  and  frequency  harmonics  of  the  problem  limit  the  numerical  so¬ 
lution  to  a  problem  such  that  with  continued  decreases  in  grid  size  and  time 
step  increment,  additional  wavelengths  and  harmonics  would  continually  be  re¬ 
vealed  while  significantly  increasing  the  computational  solution  time.  Experi¬ 
mental  techniques  in  the  measurement  of  low  velocities  and  small  temperature 
differences  resulted  in  either  velocity  measurements  being  affected  by  velocity 
sensor-induced  velocities  and  attenuated  temperature  measurements  as  a  result 
of  the  thermocouple  frequency  response.  However,  with  these  limitations  the 
numerical  and  experimental  results  obtained  demonstrate  similar  harmonics  and 
qualitative  time-dependent  behavior  of  the  developing  flow  structure. 


2.7  Conclusion* 
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Comparison  the  experimental  results  to  the  rwo-dimensional  numerical  so¬ 
lutions  revealed  the  third  dimension  of  the  laboratory  experimental  model  to 
have  a  significant  effect  upon  the  developing  flow  structure,  thus  inducing  time- 
dependent  behavior  differing  from  that  of  the  two-dimensional  numerical  solu¬ 
tions.  This  was  exhibited  in  both  the  oscillatory  behavior  at  various  chamber 
locations  as  well  as  the  decreased  amplitude  in  harmonics  of  all  the  chamber  lo¬ 
cations.  However,  there  do  appear  to  be  unique  harmonics  and  time- dependent 
behavior  at  each  of  the  chamber  locations  as  the  flow  structure  passes  through  a 
transition  of  specific  bifurcations. 

The  observed  variations  in  flow  structure  from  steady  state  to  oscillatory 
behavior  are  consistent  with  laminar  buoyant  flow  bifurcations  for  which  Prandtl 
number,  geometry,  and  Rayleig'  number  are  the  parameters  to  influence  the  flow 
structure  at  a  specific  bifurcation  (Yang,  1988).  Transition  from  a  stable,  steady 
state  flow  structure  to  one  that  is  periodic  in  time,  aperiodic,  and  finally  chaotic 
are  examples  of  specific  bifurcations  that  may  be  expected  with  variation  in 
Rayleigh  number,  Prandtl  number,  and  altered  geometry.  Such  appears  to  be  the 
case  with  the  development  of  the  flow  structure  in  the  chambers  investigated  here. 
Variation  in  Rayleigh  number  due  to  the  differences  between  the  the  internal 
heat  source  temperature,  inlet  temperature  and  wall  temperature  in  conjunction 
with  either  the  baffle  placement  or  chamber  aspect  ratios  provided  the  means  for 
these  flow  structure  transitions  to  occur.  This  resulted  in  the  observed  differences 
between  the  two-dimensional  formulation  and  laboratory  experimental  model. 
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Fig.  2.3  Two-dimensional  numerical  results  at  time  t=t„  (a)  chamber 
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Fig.  2.6  Two-dimensional  numerical  phase-space  plots  of  velocity  and 
temperature  (a)  (.501,  .31#)  (b)  (.70 1,  .31#)  (c)  (.901,  .11#)  (d)  (.501, 
.82#)  (e)  (.701,  .62#)  (f)  (.901,  .72#) 
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Fig.  2.8  Comparison  of  the  (a)  time-dependent  experimental  results  for 
velocity  and  temperature  (b)  time-dependent  numerical  results  for  velocity 
at  chamber  locations:  (a)  (.50 L,  .31//)  (b)  (.70 L,  .31//)  (c)  (.90L,  .11//)  (d) 


(.501,  .82 H)  (e)  (.70L,  .62 H)  (f)  (,90L,  .72//) 
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Fig.  2.8  Comparison  of  the  (a)  time-dependent  experimental  results  for 
velocity  and  temperature  (b)  time-dependent  numerical  results  for  velocity 
at  chamber  locations:  (a)  (.50/,  .31//)  (b)  (.70/,  .31//)  (c)  (.901,  .11//)  (d) 
(.50/,  .82//)  (e)  (.70 /,  .62 //)  (f)  (.90/,  .72//) 
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Chapter  3 


FLOW  VISUALIZATION  OF  MIXED 
CONVECTION  IN  LARGE  BAFFLED  RECTANGULAR  CHAMBERS 


3.1  Summary 


An  investigation  to  determine  the  effects  of  mixed  convection  on  the  flow 
structure  of  large  baffled  chambers  with  and  without  internal  heat  sources  has 
been  completed.  The  two-dimensional  numerical  formulation  involved  solving 
the  incompressible  time-dependent  laminar  Navier-Stokes  equations  assuming  a 
Boussinesq  fluid  with  a  Prandtl  number  of  0.7.  Flow  visualisation  in  the  cham¬ 
ber  using  2fim  Silahydrocarbcn  aerosol  droplets  generated  from  a  Collison  com¬ 
pressed  air  nebulizer  was  used  for  comparison  to  the  two-dimensional  numerical 
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solutions.  Two  baffles  were  placed  symmetrically  about  the  vertical  axis  with 
and  without  internal  heat  sources  and  asymmetric  vertical  wall  temperatures. 

Fluid  velocities  obtained  from  flow  visualization  information  using  the  Silahy- 
drocarbon  aerosol  droplets  with  and  without  internal  heat  sources  in  the  full-size 
chamber  agreed  well  with  the  predicted  velocities  from  the  two-dimensional  nu¬ 
merical  solutions.  However,  the  third  dimension  did  have  a  significant  effect  on 
the  development  of  this  flow  structure  and  the  effects  cannot  be  discounted  as 
being  negligible. 

3.2  Introduction 


Large  baffled  rectangular  chambers  with  low  inlet  velocities  have  traditionally 
been  used  to  study  a  variety  of  chemical  atmospheres  generated  for  health  effects 
and  biological  studies,  chemical  reaction  and  chemical  species  formation  stud¬ 
ies,  and  aerosol  characterization  studies.  Generally,  they  have  a  cross-sectional 
diameter  greater  than  .5  meter  with  an  inlet  velocity  set  to  maintain  10-15  cham¬ 
ber  volume  changes  per  hour.  The  inlet  and  exhaust  may  be  mounted  vertically 
or  horizontally  with  an  array  of  internal  baffles  included  for  atmosphere  disper¬ 
sion  purposes.  Additional  capabilities  include  cage  units  for  housing  laboratory 
animals  for  a  variety  biological  research  applications. 
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There  have  been  many  investigations  concerning  the  Rayleigh-Benard  flow 
in  confined  and  partially  open  enclosures  with  and  without  baffles  or  dividers 
and  internal  or  wall  heat  sources  (Bajorek  and  Lloyd,  1982;  Markatos  and  Ma- 
lin,  1982;  Powe  and  Warrington,  1983;  Nansteel  and  Greif,  1981;  Sparrow,  et 
al.,  1984;  Warrington  and  Powe,  1985;  November  and  Nansteel,  1987;  Lee  and 
Goldstein,  1988;  Ostrach  and  Austin,  1988).  Such  investigations  have  made  use 
of  both  experimental  and  numerical  results  for  information  with  regard  to  the 
temperature  profile  and  flow  structure  in  various  closed  cavity  geometries  with 
differing  magnitudes  of  Grashof,  Prandtl,  and  Rayleigh  numbers. 

Yang  (1988)  discusses  the  buoyant  enclosure  problem  as  a  dissipative  dy¬ 
namic  system  whose  flow  structure  development  is  one  of  a  series  of  bifurcations 
brought  about  by  the  driving  parameters  of  Rayleigh  number,  Prandtl  num¬ 
ber,  and  enclosure  geometry.  Development  of  a  flow  structure  from  one  that 
is  steady  state,  periodic  in  time,  aperiodic,  and  finally  chaotic  are  examples  of 
specific  bifurcations  that  may  be  expected  with  variation  in  Rayleigh  number, 
Prandtl  number,  and  altered  geometry.  Yang  (1988)  further  discusses  the  use 
of  an  attractor,  a  concept  used  in  nonlinear  dynamics,  to  describe  the  temporal 
asymptotic  behavior  of  a  trajectory  in  the  phase  space  corresponding  to  a  specific 
bifurcation. 

Initial  investigations  of  these  chambers  incorporated  flow  visualization  tech¬ 
niques  using  smoke  and  dyes  to  optimize  the  chamber  geometry  eo  that  well- 
dispersed  atmospheres  could  be  obtained  (Carpenter  and  Beethe,  1981;  Moss, 
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1981).  Quantitative  approaches  to  determine  the  mean  residence  time  of  a  cham¬ 
ber,  the  rate  of  dispersion  within  a  chamber  of  a  tracer  gas,  and  the  distribution 
within  a  chamber  using  vapors,  droplets,  and  solid  particles  have  also  been  done 
(Hemenway  et  al.,  1982;  Beethe  et  al.,  1979;  Moss,  1982;  Yeh  et  al.,  1986).  Most 
recently,  Yerkes  et  al.  (1989)  investigated  the  mixed  convection  in  large  baf¬ 
fled  rectangular  chambers  and  the  characteristic  flow  structure  development  for 
’’aided”  and  "opposed”  buoyant  forces.  Development  of  the  chamber  flow  struc¬ 
ture  was  found  to  be  predominantly  buoyant  in  nature  and  sensitive  to  small 
variations  in  the  temperature  difference  between  the  inlet  and  wall  varying  from 
steady  state  to  asymmetric  oscillatory  behavior.  Continued  investigations  ad¬ 
dressing  the  effects  of  mixed  convection  with  internal  biological  heat  sources  on 
the  on  the  transient  development  of  the  flow  structure  were  also  done  by  Yerkes 
et  al.  (1990).  However,  the  investigations  by  Yerkes  (1990)  fail  to  relate  the 
computational  numerical  results  to  the  full-size  chambers  of  interest. 

A  typical  chamber  used  for  the  purposes  described  is  shown  in  Fig.  3.1. 
The  chamber  consists  of  a  conical  inlet  intended  to  uniformly  distribute  the 
atmosphere  of  interest  across  the  inlet  plane  and  an  exhaust  manifold  to  achieve 
uniform  flow  characteristics  across  the  exhaust  plane.  Inlet  and  exhaust  flow  rates 
are  controlled  independently  to  maintain  a  slight  negative  chamber  pressure  with 
respect  to  the  ambient  pressure  (2.5  -  5.1  cm  HjO).  Biological  heat  sources  may 
consist  of  a  variety  of  laboratory  animal  species  and  are  housed  in  uniformly 
spaced  cage  units  located  directly  over  the  baffles. 
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The  scope  of  this  investigation  is  shown  in  Table  3.1  and  compares  the  two- 
dimensional  numerical  solutions  with  a  flow  visualization  technique  using  Silahy- 
drocarbon  aerosol  droplets  as  the  visual  tracer  in  the  full-size  chamber.  Of  inter¬ 
est  is  the  development  of  the  flow  structure  with  and  without  internal  biological 
heat  sources  with  asymmetric  vertical  wall  temperatures.  The  two-dimensional 
numerical  formulation  corresponds  to  a  full-size  chamber  with  a  depth  of  one 
meter.  Only  the  upper  baffles  and  heat  sources  were  considered  for  this  investi¬ 
gation. 

Biological  heat  sources  were  treated  as  a  constant  temperature  surface  at 
32°  C  and  of  a  size  and  mass  that  would  be  appropriate  for  a  200  gm  to  300  gm 
rat  in  a  confined  enclosure  (Mauderly,  1986).  It  was  assumed  that  physiolog¬ 
ical  control  of  the  body  temperature  is  sufficient  to  maintain  a  constant  body 
temperature  in  the  low  air  velocities  within  the  confines  of  the  chamber. 

3.3  Analysis 


Math  Formulation 

Two-dimensional  numerical  solutions  of  the  incompressible  time-dependent 
laminar  Navier-Stokes  equations  were  calculated  in  accordance  with  the  two- 
dimensional  problems  shown  in  Fig.  3.2a  and  3.2b  for  cases  #1  and  #2  respec¬ 
tively.  Only  the  upper  baffles  with  and  without  heat  sources  were  considered  for 
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numerical  formulation.  The  boundary  conditions  for  temperature  on  the  right 
and  left  vertical  v/alls  were  assumed  to  be  constant  and  were  obtained  from  the 
experimental  conditions.  For  this  investigation,  only  the  dimensional  form  of  the 
governing  equations  were  solved,  thereby  circumventing  the  immediate  need  for 
scaling  considerations  of  the  non-dimensional  parameters  as  is  typicallj  required 
for  buoyant  flow  conditions  (Ostrach  and  Austin,  1988). The  Boussinesq  approx¬ 
imation  was  used  in  the  vertical  direction  to  account  for  buoyancy  effects.  Using 
this  approximation  the  vertical  body  force  becomes: 


Bx  =  pg[  1  -  0(T  -  Tin)} 


where: 

for  an  ideal  gas 

For  this  investigation,  the  reference  values  for  all  of  the  gas  properties  were 
taken  to  be  those  at  the  inlet  plane. 

The  chamber  inlet  cone  was  assumed  to  provide  a  uniform  velocity  distri¬ 
bution  across  the  inlet  plane  with  the  exhaust  manifold  providing  a  velocity 
distribution  resulting  from  a  constant  exhaust  pressure  across  the  outlet  plane. 
Viscous  dissipation  and  pressure  work  were  assumed  to  be  negligible  due  to  the 
low  buoyancy-induced  velocities.  With  these  assumptions,  the  two-dimensional 
governing  equations  become: 
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Conservation  of  mass: 


dv  dw 

dy  +  ~dl  ~ 


(3.1) 


Conservation  of  momentum: 


+P90(T  -  Tin)  (3.3) 


Conservation  of  energy: 


The  boundary  conditions  were  specified  with  the  assumption  that  there  is  a 
uniform  inlet  velocity  profile  and  the  vertical  end  wall  temperature  is  constant. 
Ouly  the  magnitude  of  the  inlet  velocity  is  specified.  The  magnitude  of  the  outlet 
velocity  is  such  as  to  satisfy  the  conservation  of  mass  with  a  constant  reference 
pressure  across  the  outlet  plane. 

The  outlet  boundary  condition  for  temperature  was  considered  to  be  locally 
parabolic  such  that  the  Peclet  number  is  sufficiently  large  so  as  to  exhibit  lo¬ 
cal  one-way  behavior  in  the  axial  direction.  The  baffles  were  treated  as  being 
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infinitesimally  thick  with  a  no-slip  boundary  condition.  The  baffle  temperature 
was  specified  to  be  that  of  the  local  fluid  temperature. 


Boundary  conditions: 


Inlet:  v  =  0,  w  =  W/jy,  T  —  Tjff 


Right  wall:  v  =  0,  w  =  0,  T  =  Trwau 


Left  wall:  v  —  0,  w  =  0,  T  =  Tiwaji 


Baffle:  v  =  0,  w  —  0 


Heat  sources:  v  —  0,  w  =  0,  T  =  32°C 


Outlet:  v  =  0,  p+  =  0 


Aerosol  Science 

For  effective  use  of  flow  visualization  using  Silahydrocarbon  aerosol  droplets 
one  must  consider  some  basic  aerosol  science  to  determine  terminal  settling  ve¬ 
locities  and  the  effects  of  thermophoresis  induced  velocities  for  the  droplets  gen¬ 
erated,  (Hinds,  1982).  Droplet  behavior  in  the  flow  f  eld  will  not  necessarily  give 
correct  information  as  to  the  magnitude  of  the  fluid  velocity  due  to  the  external 
forces  acting  on  the  aerosol  droplets  resulting  in  droplet  motion  deviating  from 
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that  of  the  fluid  flow  field.  Figure  3.3  shows  the  free  body  diagram  of  a  droplet 
that  is  suspended  in  the  flow  stream.  The  vertical  forces  consist  of  the  droplet 
weight,  buoyant  force,  and  drag  force  in  which: 

F  =  0  =  -mg  +  [buoyant  force]  +  [drag  force]  (3.5) 


where 


m  =  Pdii 


[buoyant  force]  =  pfi^g 
[drag  force]  =  ZirpVd 


The  drag  force  assumes  Stokes  law  where  the  drag  coefficient,  Cd  =  24/Re. 
However,  Stokes  law  assumes  that  the  relative  velocity  of  the  fluid  at  the  surface 
of  the  droplet  is  zero.  For  droplets  less  than  1/tm  in  diameter,  the  drag  force 
will  also  include  a  slip  correction  factor,  Ce,  (Cunningham  slip  correction  factor) 
accounting  for  the  "slip”  of  a  droplet  as  the  diameter  approaches  the  mean  free 
path  of  the  fluid  molecules.  Now,  the  drag  force  becomes: 

Zirp,Vd 


[drag  force]  = 


Ce 


where 


1  + 


2.52A 


d 


(3.6) 
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for  droplet  diameters  down  to  0.1/rm. 

The  velocity  magnitude  of  the  fluid  to  maintain  the  droplet  in  a  stationary 
position  may  now  be  calculated.  This  velocity  is  the  terminal  settling  velocity  of 
a  spherical  droplet  of  a  given  density  and  diameter  where: 

v  =  <3-7) 


Thermophoresis  induced  velocity  must  alro  be  considered.  This  is  the  ve¬ 
locity  induced  by  forces  resulting  in  the  difference  in  molecular  bombardment 
on  opposite  sides  of  the  droplet  due  to  a  temperature  gradient  in  the  flow  field. 
The  magnitude  of  thermophoretic  velocity  depends  on  gas  and  droplet  proper¬ 
ties  and  on  the  temperature  gradient.  This  velocity  may  be  approximated  by  the 
following  relation  where  the  droplet  diameter  is  greater  than  the  mean  free  path 


(d  >  A): 


VtH  = 


-ZfxC'HVT 

2p/Tj 


(3.8) 


where 


(  1  \  f  kf/kd  +  iA\/i  \ 

\l  +  6A Id)  \l  +  2kf/kd  +  S.8\ld) 


(3.9) 


Motion  of  the  aerosol  droplets  may  now  be  used  and  compared  to  the  fluid 
motion.  For  instance,  if  a  droplet  has  no  vertical  motion,  the  fluid  velocity  must 
be  that  of  the  terminal  droplet  velocity  and  conversely,  if  a  droplet  is  settling  at 
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the  terminal  velocity  there  is  no  vertical  fluid  velocity.  Also,  in  regions  of  high 
temperature  gradients  such  as  at  the  vertical  wall,  the  lateral  velocity  component 
of  the  droplet  toward  the  wall  will  be  in  part  due  to  the  thermophoretic  induced 
velocity  and  in  part  due  to  droplet  diffusion. 

3.4  Numerical  Scheme 


The  numerical  scheme  used  is  a  finite-difference  iterative  method  of  solution 
using  a  control-volume  approach  as  developed  by  Spalding  et  al.  (1880).  A 
two-dimensional  numerical  model  was  formulated  using  the  fully  implicit  finite- 
difference  scheme  of  the  time  dependent,  laminar  Navier-Stokes  equations.  The 
"SIMPLEST”  (Spalding  et  al.,  1980)  method  of  solution  for  the  momentum 
equations  was  used  with  the  hybrid  differencing  formulation.  When  a  cell  Peclet 
number  is  within  the  range  -2  to  2,  a  central-difference  scheme  is  used  and  when 
the  cell  Peclet  number  is  outside  this  range  the  upwind  differencing  scheme  is 
used. 

The  outlet  boundary  condition  for  temperature  was  considered  to  be  locally 
parabolic  such  that  the  Peclet  number  is  sufficiently  large  so  as  to  exhibit  local 
one-way  behavior  in  the  axial  direction  (Patankar,  1980).  This  approach  to 
specifying  the  outlet  boundary  condition  for  a  dependent  variable  is  generally 
accepted  as  the  most  available  approach  for  internal  flow  problems  in  which 
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the  fluid  leaves  the  calculation  domain  and  the  dependent  variable  of  interest  is 
unknown. 

The  baffles  were  treated  as  being  infinitesimally  thick  with  a  no-slip  bound¬ 
ary  condition  by  setting  the  cell  wall  velocity  to  zero  o'rer  the  region  of  the  baffle 
location.  The  baffle  temperature  was  specified  to  be  that  of  the  local  fluid  tem¬ 
perature.  Heat  sources  were  ’’blocked”  by  setting  the  cell  velocity  components 
within  the  heat  source  region  to  zero  and  fixing  to  a  constant  temperature. 

The  solution  sequence  involved  solving  for  the  velocity  and  temperature  fields 
from  an  assumed  pressure  field  at  a  specific  time  step.  The  pressure  field  was 
then  subsequently  updated  using  these  velocity  and  temperature  fields  in  the 
pressure-correction  equation  such  that  the  conservation  of  mass  was  satisfied. 
This  iterative  sweep  process  was  then  repeated  using  a  slab-by-slab  method  at 
each  time  step  until  convergence  criteria  were  satisfied.  Convergence  charac¬ 
teristics  at  each  time  step  were  mai;  »ained  by  veducing  the  time  step  where 
appropriate  with  the  minimum  computational  time  required  to  reach  a  steady 
state  solution  being  at  least  three  characteristic  times,  H/Wjff,  or  approximately 
ten  minutes.  In  addition,  conditions  for  convergence  to  a  solution,  either  at  a 
specific  time  step  or  to  a  steady  state  solution,  was  based  upon  a  dependent 
variable,  <f> ,  varying  less  than  a  predetermined  change  of  the  magnitude  fraction, 
e^,  between  successive  sweeps  or  time  steps  where: 


I  <£«  -  <P»+ 1 !  !4>  < 
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The  values  of  were  varied  and  it  was  required  that  the  magnitude  of  the 
dependent  variable  be  greater  than  the  roundoff  error  or  truncation  error  such 
as  would  happen  when  velocities  approach  zero.  Typically,  for  a  steady  state 
solution,  ranged  from  1.0  x  10“ 3  to  1.0  x  10“ 2  over  a  time  span  of  25  seconds, 
while  ever  five  successive  sweeps  was  less  than  1.0  x  10“  3 .  Intermediate  results 
were  monitored  to  insure  convergence  and  numerical  stability  at  a  specific  time 
step  during  the  course  of  the  solution.  These  intermediate  results  were  also 
used  to  determine  convergence  toward  either  a  steady  state  or  a  steady  periodic 
solution. 

Appropriate  grid  size  and  spacing  were  cl  osen  to  give  acceptable  numerical 
accuracy  while  still  maintaining  reasonable  and  acceptable  computational  times. 
Due  to  the  complexity  of  the  flow  structure,  the  errors  as  a  result  of  grid  size  were 
determined  first,  by  comparing  the  shape  of  the  flow  structure  for  similarity  with 
a  previous  solution  at  a  larger  grid  size.  Secondly,  random  locations  throughout 
the  chamber  were  selected  and  the  average  error  compared  between  the  similar 
solutions  such  that  decreasing  the  grid  size  by  ten  percent  resulted  in  less  than  a 
one  percent  error  in  the  dependent  variable.  For  this  problem,  a  uniform  100  x  100 
grid  was  chosen  to  give  good  coverage. 
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3.5  Experimental  Approach 


Figure  3.1  shows  tne  full-size  whole-body  inhalation  exposure  chamber  as 
is  used  in  the  toxicity  testing  of  various  chemical  compounds.  The  rectangular 
chamber  body  measures  0.813  m  x  0.813  m  x  0.686  m  and  is  constructed  of 
stainless  steel  with  the  front  and  sides  fabricated  of  glass  for  observation  during 
chamber  operation.  The  inlet  consists  of  &  annulus  and  conical  section  for  the 
purpose  of  obtaining  a  uniform  velocity  distribution  across  the  inlet  plane.  The 
exhaust  consists  of  a  radial  spider  also  for  the  purpose  of  obtaining  uniform 
exhaust  characteristics  across  the  exhaust  plane. 

Figure  3.4  shows  the  schematic  of  the  aerosol  generation  apparatus  that  was 
used  to  generate  a  Silahydrocarbon  aerosol  droplet  for  the  purposes  of  visualiz¬ 
ing  the  chamber  flow  structure.  The  generation  apparatus  consisted  of  a  six-hole 
Collison  (BGI,  Inc.  Waltham,  MA.)  compressed  air  nebulizer  operated  at  20  psi 
to  produce  a  mass  median  droplet  diameter  of  2.0 /im  with  a  geometric  standard 
deviation  of  2.0.  The  droplet  density  was  assumed  to  be  thal  of  the  bulk  fluid 
density  of  Silahydrocarbon,  pj=.3145  gm/cc,  resulting  ;n  a  terminal  droplet  ve¬ 
locity  of  9.85  x  10~J  cm/s.  The  thermophoretic  velocity  for  a  1.0  fim.  droplet 
is  on  the  order  of  1.0  x  10-<  cm/s  for  a  one  degree  temperature  per  centimeter 
gradient. 
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Chamber  inlet  and  wall  temperatures  were  obtained  using  iron-constantan 
surface  mounted  thermocouples  placed  along  the  centerline  of  the  vertical  end 
walls  at  .15 H  and  .70 H.  These  thermocouples  were  calibrated  to  within  ±0.1°C 
using  a  constant  temperature  stirred  water  bath  and  an  NBS  traceable  thermome¬ 
ter  ±0.05°C.  The  vertical  wall  temperatures  were  averaged  to  give  the  temper¬ 
ature  boundary  condition  that  was  subsequently  used  in  the  two-dimensional 
numerical  formulation. 

Flow  visualization  was  obtained  using  an  8  mm  video  camera  (Sony  CCD- 
V9)  with  lighting  provided  by  insulated  florescent  lamps  placed  at  the  top  and 
base  of  the  right  vertical  end  wall.  The  fluid  velocity  in  the  plane  of  the  camera 
was  calculated  by  observation  of  the  aerosol  path  with  known  distance  references 
in  the  chamber.  Effects  of  the  lamps  were  determined  prior  to  the  experimental 
procedure  by  visual  observation  using  low  level  lighting  and  found  to  be  negli¬ 
gible.  The  experimental  procedure  consisted  of  introducing  a  30  second  bolus 
of  Silahydrocarbon  droplet  aerosol  while  maintaining  the  volumetric  flow  rate 
through  the  chamber  to  produce  an  inlet  velocity  of  0.71  x  i0-3  m/s. 

3.0  Results  and  Discussion 


Experimental  flow  visualization  information  was  obtained  for  cases  #1  and 
#2  with  the  right  and  left  vertical  end  wall  temperatures  incorporated  as  bound¬ 
ary  conditions  for  the  subsequent  corresponding  numerical  formulations  and  so- 
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lutions  for  the  cases  investigated,  Table  3.1.  Figures  3.5  and  3.6  show  the 
two-dimensional  numerical  results  for  the  velocity  direction  and  magnitude  pro¬ 
files  for  cases  #1  and  #2  respectively. 

Also  shown  in  Figs.  3.5  and  3.6  are  the  chamber  locations  for  which  the  veloc¬ 
ity  obtained  from  the  visualization  experiments  are  compared  to  the  numerical 
results.  The  chamber  locations  for  the  velocity  comparison  in  case  #1  were 
taken  at  the  lower  vertical  surface  at  the  left  wall,  A,  the  right  lower  horizontal 
plane  just  above  the  exhaust  manifold,  B,  and  the  lower  vertical  surface  at  the 
right  wall,  C.  The  chamber  location  for  the  velocity  comparison  in  case  #2  were 
taken  at  the  location  between  the  two  baffles  and  heat  sources  along  the  cham¬ 
ber  centerline,  A,  and  the  right  lower  horizontal  plane  just  above  the  exhaust 
manifold,  B.  Case  #2  also  showed  a  separation  point  located  at  the  lower  right 
vertical  wall,  C.  The  vertical  location  of  this  separation  point  is  compared  for 
both  the  numerical  results  and  measurements  made  from  the  flow  visualization 
information. 

A  summary  of  the  two-dimensional  numerical  results  compared  with  the  ex¬ 
perimental  flow  visualization  are  shown  in  Table  3.2.  The  numerical  solutions 
showed  steady  state  behavior  for  the  two  cases  studied.  There  was  good  agree¬ 
ment  between  the  numerical  and  experimental  results  for  the  chamber  locations 
investigated.  In  general,  even  though  the  full-size  chamber  geometry  differs  from 
that  of  the  numerical  formulation,  the  visual  observations  indicated  numerous 
similarities  in  the  development  of  the  flow  structure. 
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Case  #1  experimental  observations  showed  a  rotational  mass  flux  about  the 
perimeter  of  the  chamber  as  did  the  numerical  solution.  There  was  little  vertical 
motion  in  the  central  lower  region  beneath  the  baffles.  Comparison  of  the  nu¬ 
merically  obtained  vertical  velocities  in  this  lower  region  indicated  the  velocity 
magnitude  to  be  approximately  that  of  the  droplet  settling  velocity.  The  verti¬ 
cal  velocity  gradually  increased  in  magnitude  as  the  right  vertical  end  wall  was 
approached. 

Benard  cells  were  observed  for  case  #2  above  the  upper  baffles  and  heat 
sources  while  the  lower  region  below  the  baffles  showed  stagnate  layers  of  aerosol 
droplets.  Numerically  obtained  vertical  velocities  in  this  lower  region  also  showed 
the  velocity  to  be  on  the  order  of  the  droplet  settling  velocity.  Separation  of  the 
flow  from  the  lower  right  vertical  end  wall  was  due  to  a  mass  flowing  down  the 
right  vertical  end  wall  intersecting  with  mass  flowing  across  the  exhaust  plane 
and  up  the  same  vertical  end  wall.  The  intersection  location  resulted  in  the 
observed  separation  of  the  mass  from  the  vertical  wall.  The  bulk  of  the  aerosol 
mass  remained  in  the  upper  region  showing  a  significant  increase  in  residence 
time  directly  influenced  by  the  heat  sources.  An  apparent  distinguishing  feature 
of  this  case  is  the  lack  of  flow  between  the  baffles  along  the  chamber  centerline. 
The  absence  of  any  vertical  mass  transport  in  this  region  was  confirmed  by  both 
numerical  and  the  experimental  visualization  information. 

The  flow  structure  for  case  #1  appeared  to  be  sensitive  to  external  temper¬ 
ature  influences  as  was  indicated  when  a  heat  gun  was  used  to  heat  the  third 
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dimension  front  vertical  surface.  The  flow  structure  for  c ase  #2  did  not  appear 
to  be  as  sensitive  to  similar  external  temperature  influences. 

3.7  Conclusion 


The  two-dimensional  numerical  solutions  agreed  well  with  the  visual  stud¬ 
ies  that  were  conducted  for  the  cases  with  and  without  internal  heat  sources 
and  asymmetric  vertical  end  wall  temperatures.  There  were  significant  third 
dimensional  effects  that  need  to  be  considered  for  future  studies,  however,  the 
two-dimensional  solutions  gave  accurate  information  with  regard  to  the  general 
form  of  the  flow  structure.  Characteristic  elements  of  the  flow  structures  for 
both  cases  were  shown  to  be  present  in  the  numerical  solutions  as  well  as  the 
experimental  observations  for  the  chamber  locations  investigated.  Differences  in 
the  chamber  geometry  and  the  numerical  formulation  for  the  most  part  did  not 
appear  to  greatly  effect  the  results  except  for  the  inlet  cone  region  where  the  flow 
was  shown  to  follow  the  inlet  cone  surface.  The  effects  of  internal  heat  sources 
were  consistent  and  showed  a  significant  deviation  in  the  flow  structure  from  one 
without  the  internal  heat  sources.  This  deviation  in  the  flow  structure  appeared 
to  result  not  only  in  the  development  of  Benard  cells  but  also  resulting  in  the 
noticeable  increase  in  mass  residence  time  in  the  upper  region  above  the  baffles 
and  heat  sources.  For  both  cases  mass  was  transported  along  the  chamber  walls 
verifying  the  existence  of  a  high  velocity  boundary  layer. 
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Table  3.1 

EXPERIMENTAL  CONDITIONS  FOR 
2-D  NUMERICAL  COMPARISON 


case  # 

1 

W/O  HEAT  SOURCES 


2 

W  HEAT  SOURCES 


Temperature  specification 

T/jv=21.5°C 

TrWi=22.4°C 

TlmM=2l.2°C 

T/jv=24.0°C 
Trw*U= 23.0°C 
Tlw4a=22.0°C 


(upper  baffles  and  heat  sources) 


WIN  =  0.71  x  lO-’m/s 
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Table  3.2 

EXPERIMENTAL  FLOW  VISUALIZATION 
AND  2-D  NUMERICAL  COMPARISON 


case  # 

Location 

Numerical 

Experimental 

1 

A 

2-5  cm/s-IJ. 

2-3  cm/slj- 

W/O  HEAT  SOURCES 

B 

1-3  cm/s=> 

4.5  x  10"J  -2.0  x  lO^cm/sft 

4-5  cm/s=> 

C 

1-3  cm/sft 

1.0-1. 5  cm/sff 

2 

A 

9.2  cm/s=>  (max) 

10-12  cm/s=> 

W  HEATSOURCES 

B 

1.5-3. 0  cm/»=> 

5.0  x  10-*  -  2.5  x  10-acm/sft 

1-2  cm/s=> 

C 

8. 2-9. 6  cml 

5-10  cmf 

EXHAUST 
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Fig.  3.3  Force*  on  aerosol  droplet 


COMPRESSED 


OVERALL  SUMMARY 


Conclusions 


Development  of  the  chamber  flow  structure  with  and  without  internal  heat 
sources  was  shown  to  be  predominately  buoyant  in  nature  with  the  velocity  pro¬ 
file  developing  primarily  due  to  temperature  differences  within  the  confines  of  the 
chamber.  The  chamber  flow  structure  without  internal  heat  sources  was  found 
to  be  sensitive  to  small  variations  in  the  temperature  difference  between  the  inlet 
and  vertical  wall  with  an  increase  in  the  mass  flux  to  the  vertical  wall  for  inlet 
temperatures  greater  than  the  vertical  wall  temperatures.  As  the  inlet  to  vertical 
wall  temperatures  increased,  there  was  shown  ‘  be  an  increasing  low  velocity 
region  depicted  by  the  stratified  temperature  gradient  for  which  mass  and  energy 
is  transported  by  diffusion  as  opposed  to  convection.  When  inlet  temperatures 
became  less  than  the  vertical  wall  temperature  and  without  internal  heat  sources, 
the  flow  structure  became  oscillatory  in  nature  primarily  due  to  the  placement 
of  the  lower  baffles.  This  oscillatory  flow  structure  consisted  of  a  central  core 
of  cooler  fluid  falling  and  subsequently  rising  along  the  warmer  vertical  walls. 
With  the  addition  of  internal  heat  sources,  the  flow  structure  developed  multiple 
Benard  cells  typically  common  to  buoyant  flow  structures.  Development  of  the 
flow  structure  was  shown  to  be  marginally  stable  with  self- sustained  oscillatory 
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behavior  as  the  Benard  cells  formed  and  subsequently  collapsed.  The  flow  struc¬ 
ture  developed  through  a  series  of  bifurcations  from  steady-state  to  periodic, 
aperiodic,  and  chaotic  with  increasing  temperature  differences  between  the  heat 
sources,  inlet,  and  vertical  walls. 

Experimental  efforts  using  both  the  experimental  laboratory  model  and  the 
large-*,  full-size  chamber  provided  useful  information  in  comparing  the  exper¬ 
imental  results  to  the  numerical  solutions.  Use  of  the  laboratory  experimental 
model  to  depict  the  larger,  full-size  chamber  showed  significart  three-dimensional 
effects  influencing  the  development  of  the  flow  structure  within  the  model.  How¬ 
ever,  even  with  these  differences,  formulation  of  a  three-dimensional  numerical 
approach  successfully  addressed  these  inconsistencies.  Flow  visualization  in  the 
larger,  full-size  chamber  agreed  well  with  the  two-dimensional  numerical  formu¬ 
lation.  The  fluid  velocity  as  determined  at  the  selected  chamber  locations  agreed 
well  with  the  numerical  results  as  did  the  location  of  the  separation  point  when 
the  internal  heat  sources  were  considered.  There  still  appeared  to  be  significar  t 
third-dimensional  effects  in  the  full-size  chamber. 

The  two-dimensional  formulation  for  either  the  laboratory  model  or  larger, 
full-size  chamber  appeared  to  be  limited  due  to  the  neglecting  of  the  third- 
dimension.  However,  the  two  dimensional  results  did  predict  the  general  cham¬ 
ber  behavior  as  shown  by  the  visualization  experiments  in  the  full-size  chamber. 
The  assumption  of  constant  and  uniform  wall  temperatures,  heat  source  tempera¬ 
tures,  and  inlet  velocity  appeared  to  be  appropriate  for  these  initial  investigations 
as  was  also  indicated  by  the  visualization  experiment. 
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Recommendations 


With  these  initial  studies  concluded,  it  is  appropriate  to  back  track  and  ini¬ 
tiate  detailed  parametric  studies  to  determine  the  effects  of  temperature,  baffle 
placement  or  altered  geometry,  and  heat  source  placement.  Coinciding  with  these 
parametric  studies,  evaluation  of  the  multiple  length  scales  for  future  investiga¬ 
tions  would  be  helpful  in  addressing  the  order  of  magnitude  effects  of  various 
non-dimensional  terms  at  various  chamber  locations  in  addition  to  deriving  the 
non-dimensional  and  normalized  sets  of  governing  equations. 

These  chambers  are  normally  operated  with  multiple  chemical  species,  mul¬ 
tiple  phase  atmospheres  with  and  without  condensation,  chemical  reactions,  etc. 
Therefore  additional  efforts  in  investigating  the  transport  characteristics  of  a  va¬ 
riety  of  airborne  compounds  is  of  interest  in  further  developing  and  describing 
the  physics  of  chamber  performance  under  specific  applications. 

Although  the  computational  time  for  the  solution  of  these  problems  has  been 
extensive,  an  expanded  formulation  to  one  that  incorporates  the  effects  of  three- 
dimensions  is  desirable.  These  formulations  would  more  accurately  depict  cham¬ 
ber  behavior  defining  criteria  such  as  inlet  and  exhaust  configurations  for  future 
chamber  design  and  modification. 
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Appendix  A:  BUOYANT  FLOW  AND  EXPOSURE  CHAMBER  OVERVIEW 

A)  BUOYANT  FLOW  OVERVIEW 

These  investigations  involve  the  development  of  a  velocity  field  within  a  fluid 
primarily  due  to  changes  in  the  fluid  density  as  a  result  of  temperature  gradients. 

It  is  appropriate  to  discuss  some  basic  theory  behind  the  development  of  this  type 
of  fluid  motion  to  further  supplement  the  investigations  presented  herein. 

It  may  be  recalled  that  mass  and  energy  is  transported  either  by  convection 
(velocity)  or  diffusion  (random  molecular  motion).  Convection  transport  ci  n 
be  induced  by  either  external  means  such  as  fluid  machinery  or  natural  mea  is 
such  as  gravitational,  centrifugal,  or  coriolis  forces.  Motion  is  induced  when  he 
lighter  less  dense  fluid  rises  and  the  heavier  greater  density  fluid  falls. 

There  have  been  numerous  investigations  into  the  development  of  Rayleigh- 
Benard  flow  within  confined  enclosures.  However  there  is  little  information  con¬ 
cerning  the  addition  of  low  velocity  fluid  flow  into  these  same  enclosures.  Figure 
A.i  shows  the  scope  of  previous  investigations  for  free,  forced,  and  mixed  con¬ 
vection  regimes  for  flow  in  vertical  tubes,  (Aung,  1987).  There  has  been  very 
little  work  with  low  Reynolds  number  and  high  Rayleigh  number  flows  in  cavities 
whose  diameter  to  height  ratio  approach  one.  For  the  investigations  conducted 
herein,  inlet  Reynolds  numbers  range  from  32  to  235  with  Rayleigh  numbers  on 
the  order  of  107.  Although  the  chamber  of  interest  is  of  a  rectangular  cross- 
section  rather  than  a  circular  cross-section,  the  hydraulic  diameter  to  chamber 


\ 
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height  ranged  from  0.18  to  1.31  for  the  experimental  laboratory  model  and  full- 
size  chamber  respectively.  The  scope  of  this  investigation  therefore  encompasses 
the  free  convection,  laminar  flow  regime  for  which  there  has  been  very  little 
investigative  research  efforts. 

The  general  form  of  the  time- dependent  differential  governing  equations  for 
a  compressible  variable  property  fluid  are  as  follows: 

Conservation  of  mass: 


dp 

at 


+  v .  (Pv)  =  o 


Conservation  of  momentum: 


(A.1) 


DV 
P  Dt 


=  V-r 


Vp  +  pB 


Conservation  of  energy: 


(A.2) 


„C,— =  V-Jt  +  S  M.3) 

where  D/Dt  is  the  substantia!  derivative  (d/dt  +  V*  •  V).  The  energy  equation 
(.4.3)  does  not  contain  the  viscous  dissipation  and  pressure  work  terms,  since  both 
are  small  due  to  the  low  velocities  encountered  in  natural  convection  phenomena 
(Gebhart,  1971). 
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Pa 


For  these  investigations,  the  fluid  property  variation  vas  accounted  for  by 
assuming  the  fluid  to  be  a  Boussinesq  fluid.  This  approximation  relates  the  fluid 
density  in  the  body  force  term  to  a  reference  fluid  density  by  way  of  the  fluid 
expansion  coefficient,  /3,  and  difference  between  the  fluid  temperature  and  a  ref¬ 
erence  fluid  temperature  (Yang,  1987;  Gebhart,  1971).  By  Bernoulli’s  equation, 
the  local  depression  of  pressure  in  the  flow  region  at  any  given  vertical  distance, 
z,  is  not  greater  than  pu;2/ 2  or 

Po~P<  0[gz(p0  -  p)\ 


Using  this  relation  and  since  p  =  p(p,t ),  the  differences  in  density  at  a  given 
z  may  be  written  as  a  double  expansion  in  terms  of  the  temperature  difference 
and  pressure  difference  where: 


Using  the  estimate  for  pa  -  p,  the  density  difference  for  an  ideal  gas  (p  =  pRT), 
in  terms  of  (T  —  T0)  becomes 


P0{T  -  T0)  90*.  u#/,/T  Tw  N 

Po  P  1  +  (3{T  -  T0)  +  R^Po  P^  +  R  ^  °^P°  P> 


-(-additional  cross  terms 


where 
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for  an  ideal  gas 


The  quantity  g(3z/ R  is  very  small  for  gases  at  the  terrestrial  intensity  of 
gravity  and  (T  —  T0)/T  is  generally  assumed  to  be  small.  Now,  the  following 
approximation  for  density  becomes 

p  =  po[l-0{T-To)]  ( AA ) 

Using  this  approximation  for  the  variation  of  density  in  the  body  force  term 
of  the  momentum  equations,  the  formulation  becomes  one  in  terms  of  the  differ¬ 
ence  between  the  fluid  temperature  and  a  reference  fluid  temperature.  Allowing 
all  the  remaining  density  terms  in  the  governing  equations  to  be  the  constant 
reference  fluid  density,  p  =  p0,  the  governing  equations  can  now  be  formulated  in 
the  incompressible  form.  As  with  any  approximation,  there  are  limits  for  which  a 
Boussinesq  fluid  may  be  assumed  as  discussed  by  Zhong  (1985).  This  approxima¬ 
tion  is  generally  valid  for  9a  <  0.1  with  an  overprediction  of  the  vertical  velocity 
as  much  as  20  percent  as  9„  approaches  0.2.  For  this  investigation,  the  maximum 
value  of  90  is  .027  which  is  well  within  the  limits  of  acceptable  use  for  assuming 
a  Boussinesq  fluid. 

Using  this  approximation  in  the  vertical  direction  to  account  for  buoyancy 
effects,  the  vertical  body  force  becomes: 

Bt=pg  =  pog[l-I3(T-T0))  (A.5) 


/ 
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The  p0g  term  in  the  Boussinesq  approximation,  equation  (j4.2),  may  now  be 
transformed  into  a  hydrostatic  pressure  term  an d  be  combined  with  the  existing 
static  pressure  term  in  the  momentum  equation  while  noting  the  gravitational 
direction  to  maintain  the  correct  sign  in  the  body  force  term.  For  these  investi¬ 
gations,  the  reference  values  for  all  of  the  gas  properties  were  taken  to  be  those 
at  the  inlet  plane.  Viscous  dissipation  and  pressure  work  were  assumed  to  be 
negligible  due  to  the  low  buoyancy-induced  velocities.  With  these  assumptions, 
the  time-dependent  incompressible  governing  equations  become: 

Conservation  of  mass: 


(A.  6) 


Conservation  of  momentum: 


ac-direction: 

'  (§? +  v '  =  +  M-7) 


y-direction: 

',(^  +  v-H)  =  -^-+v'(',v',)  (as) 


2-direction: 


V 

di 


+  V-(MVu ,)  +  pS0{T-T.) 


(A.9) 
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Conservation  of  energy: 


'(f +v- (<*))-"  (£vr)  <A10> 

Additional  difficulties  arise  that  are  specific  to  buoyant  flow  problems  in  con¬ 
fined  enclosures  dealing  with  the  multiple  length  and  velocity  scales  (Ostrach 
and  Austin,  1988).  Non-dimensionalizing  and  normalizing  the  dimensional  form 
of  the  governing  equations  require  physically  correct  reference  lengths  and  char¬ 
acteristic  velocities  to  accurately  solve  the  problem  while  still  maintaining  the 
physics  of  the  problem.  With  this  specific  problem,  as  with  other  buoyant  flow 
problems,  there  are  multiple  scale  lengths  and  characteristic  velocities  depending 
upon  the  chamber  location.  The  physics  defining  the  fluid  mechanics  differ  due 
to  the  prevailing  dominate  forces  at  various  chamber  locations.  These  dominate 
forces  may  vary  from  that  of  buoyancy  forces,  viscous  forces,  and  diffusion  forces. 
To  accurately  solve  and  reflect  the  physics  of  the  entire  problem,  the  different 
reference  length  scales  and  characteristic  velocities  in  each  of  these  regions  must 
be  correctly  identified  and  be  reflected  in  the  non-dimensional  and  normalized 
sets  of  governing  equations.  Since  these  first  sets  of  investigations  were  for  the 
purpose  of  identifying  the  physics  of  the  problem,  only  the  dimensional  sets  of 
governing  equations  were  solved  numerically.  With  these  initial  and  further  para¬ 
metric  studies,  the  appropriate  development  of  the  problem  in  a  non-dimensional 
normalized  form  may  be  addressed  while  still  accurately  addressing  the  physics 
of  the  problem. 
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This  difficulty  can  best  be  illustrated  by  considering  the  possible  multiple 
characteristic  velocities  in  terms  of  the  characteristic  length  scale  depending  upon 
the  dominate  forces  governing  the  physics  of  the  problem,  (Ostrach,  1982): 

Vc  =  v/Lc  =  (viscous  vs.  inertia) 

Ve  =  a/ Lc  (convection  vs.  conduction) 

Vc  =  (av)  s/Le  (combination) 

Vc  =  g0&TcL\]v  (buoyancy  vs.  viscous) 


Ve  =  ( g(3ATeLc)' 5  (buoyancy  vs.  inertia) 
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B)  EXPOSURE  CHAMBER  OVERVIEW 

Research  efforts  to  determine  the  inhalation  toxicity  of  a  variety  of  airborne 
compounds  utilize  either  whole-body  or  nose-only  exposure  methods.  As  a  re¬ 
sult,  airborne  compounds  or  toxicants  that  may  consist  of  vapour,  droplets,  solid 
particulates,  or  a  combination  of  the  three  have  to  be  transported  with  maximum 
efficiency  through  a  variety  of  flow  elements  so  as  to  be  inhaled  by  experimental 
animals.  Transport  tubes,  mixers,  and  exposure  chambers,  for  example,  have  in¬ 
herent  loses  associated  with  their  operation  during  the  course  of  an  experiment. 
This  study  addresses  the  whole-body  exposure  chamber  and  uses  computational 
methods  in  conjunction  with  experimental  data  to  develop  a  better  understanding 
of  physics  of  chamber  operation  for  improvement,  modification,  and  future  de¬ 
velopment.  Transport  phenomenon  such  as  diffusion  and  the  convection  of  mass, 
enthalpy,  and  momentum  play  an  important  role  on  the  efficiency  of  chamber 
operation.  With  the  addition  of  heat  generation  from  animal  species  within  the 
chamber,  these  transport  phenomenon  are  further  complicated  with  the  intro¬ 
duction  of  the  free  convection  of  heat,  possibly  altering  flow  patterns  within  the 
chamber. 

There  have  been  several  exposure  chamber  configurations  developed  for  in¬ 
halation  toxicology  studies.  Figure  A. 2  shows  four  examples  of  the  variety  of 
possible  chamber  configurations  which  have  been  or  are  being  used  for  inhalation 
research.  Figure  A. 2a  shows  a  horizontal  flow  chamber  developed  by  Hemenway 
(1982)  for  which  the  toxicant  is  dispersed  horizontally.  Figure  A.2b  shows  the 
multi-tiered  chamber  developed  by  Moss  (1981)  while  Fig.  A.2c  and  A.2d  show  a 
conventional  pyramidal  end  cone  "Hinners-type"  (Carpenter  and  Beethe,  1981) 
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and  cylindrical  chamber  (Riley,  1986)  type  configurations.  These  chambers  axe 
generally  fabricated  of  stainless  steel  and  glass  in  order  to  minimize  toxicant  at¬ 
mospheres  reacting  with  the  chamber.  The  primary  purpose  of  these  chambers 
has  been  to  deliver  will  dispersed  toxicants  at  low  velocities  to  determine  the 
toxic  effects  of  various  compounds. 

Evaluation  of  these  chambers  have  been  performed  by  numerous  investigators 
addressing  spatial  and  temporal  toxicant  distributions  within  *he  chambers  of 
interest  for  specific  applications  (Beethe  et  al.,  1979;  Yeh  et  al.,  1986).  The 
approach  to  addressing  these  spatial  and  temporal  variations  has  been  to  sample 
at  predetermined  chamber  locations  by  extracting  a  volume  of  air  and  toxicant 
over  a  given  period  of  time  and  subsequently  comparing  the  variation  of  the 
toxicant  distribution  within  the  chamber. 

Additional  evaluations  of  these  chambers  using  flow  visualization  and  quanti¬ 
tative  models  addressing  the  mean  residence  time  of  the  chamber  and  the  rate  of 
dispersion  of  a  toxicant  within  the  chamber  have  also  been  performed  (Carpen¬ 
ter  and  Beethe,  1981;  Moss,  1981;  Hemenway  et  al.,  1986).  Flow  visualization 
techniques  by  Moss  (1981)  utilized  a  hydraulic  water  model  with  the  injection  of 
a  tracer  dye  as  the  visual  reference  in  determining  the  flow  structure  within  the 
chamber  while  Carpenter  (1981)  injected  TiOj-HCL  smoke  into  the  chamber  to 
ascertain  the  chamber  flow  structure.  Quantitative  information  for  the  chamber 
as  a  whole  was  developed  by  Hemenway,  et  al.  (1986)  by  addressing  the  chamber 
response  to  an  injection  of  a  toxicant  or  test  material  by  comparing  the  inlet  and 
exhaust  profiles  over  time.  This  gave  information  as  to  the  degree  of  chamber 
mixing,  dead  volume  size,  and  shunting  of  the  test  material  to  the  exhaust. 
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The  effects  of  internal  heat  sources  or  animal  species  have  been  addressed  by 
Bernstein  and  Drew  (1981)  and  Yeh  et  al.  (1988).  Bernstein  and  Drew  (1981) 
discussed  animal  loading  and  the  effectiveness  of  eliminating  or  dissipating  the 
heat  load  from  the  chamber.  The  effects  of  the  heat  sources  on  the  chamber  as 
a  whole  is  discussed  by  Yeh  (1988)  showing  the  effects  of  spatial  an  temporal 
variation  with  and  without  animals. 

Although  these  investigations  have  provided  useful  insight  and  information 
as  to  the  overall  performance,  to  date,  there  has  been  little  use  made  of  com¬ 
putational  fluid  dynamics  for  the  evaluation  and  design  of  these  chambers,  nor 
have  consideration  as  to  what  role  the  buoyancy  forces  play  and  what,  if  any 
are  the  possible  effects  due  to  sampling  in  determining  spatial  and  temporal 
distributions. 
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Figure  A.l  Free,  forced,  and  mixed  convection  regimes  for  flow  in  vertical 
circular  tubes  for  uniform  heat  flux  (UHF)  and  uniform  wall  temperature 


Fig.  A. 2  Posrble  chamber  configurations  in  use  for  inhalation  research 
(a)  horizontal  flow  chamber  (b)  multi-tiered  chamber  (c)  pyramidal  end  cone 
"Hinners-type”  chamber  (d)  cylindrical  chamber 


Appendix  B:  EXPERIMENTAL  INSTRUMENTATION 


A)  Rotameter  pressure  correction,  calibration  and  use 

i)  Theory  and  use 

The  volumetric  air  flow  rate  through  the  experimental  chamber  was  moni¬ 
tored  and  controlled  using  a  Matheson  "tube  cube”  rotameter.  It  is  appropriate 
to  review  the  theory  of  operation  for  these  devices  in  addition  to  a  method  for 
pressure  correction  to  determine  the  actual  volumetric  flow  rate  during  operation. 

A  rotameter  consists  of  a  tapered  tube  in  which  some  sort  of  a  "float”  assumes 
a  vertical  position  corresponding  to  a  given  flow  rate  through  the  tube.  The 
tapered  tube  diameter  varies  continuously  with  vertical  displacement  and  serves 
three  purposes: 

1)  Balance  of  the  float  is  self  maintaining. 

2)  Provides  a  constant  presaure  drop  ( 8P )  for  different  flow  rates.  [Since 
6 P  =  6P(Q J),  for  a  fixed  area,  the  tapered  tube  provides  a  variable  area 
such  that  the  6P  is  constant.] 

3)  If  the  taper  varies  linearly  with  height,  the  float  position  will  vary  linearly 
with  flow  rate. 
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Figure  B.l  shows  the  free  body  diagram  of  a  ball  float  that  is  suspended 
in  the  flow  stream.  The  vertical  forces  consist  of  the  ball  weight,  buoyant  force, 
and  drag  force  in  which: 

^2  F  —  0  =  ~mg  +  [buoyant  force]  +  [drag  force]  (J3.1) 

where: 


m  -  PfcV), 

[buoyant  force]  = 
idr^s  force]  =  V 

Solving  equation  (B.l)  for  velocity  and  substituting  into  the  conservation  of 
mass  where: 


(£.2) 


The  cross-sectional  area,  Ay  can  further  be  deflned  in  terms  of  tube  taper 
and  vertical  position,  y  (note:  D0  =  d): 

A  =  -((£>»  +  ay)1  -  d* )  as  ^D9ay  »  ^da y 


(B.  3) 
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Substituting  equation  (5.3)  into  (5.2)  and  letting 


K  = 


(  2V6g  \ 

U  AhCD)J 


5 


we  get: 

=  Ky[(pb  -  Pf)pf}  5  {ba) 

Now  we  have  a  relationship  between  the  mass  flow  rate  and  the  vertical 
displacement  of  the  ball.  For  viscous  fluids,  the  drag  coefficient  can  be  made 
nearly  independent  of  viscosity  by  using  sharp  edge  floats.  For  an  ideal  gas, 
in  which  a  change  in  density  will  occur  at  differing  pressures  by  way  of  the 
ideal  gas  law,  there  must  be  a  pressure  correction  made  to  maintain  an  accurate 
measurement  of  the  volumetric  flow  rate  while  the  rotameter  is  in  use. 


ii)  Rotameter  pressure  correction 

Consider  two  application  in  which  the  first  is  the  rotameter  operating  at  some 
reference  pressure  and  temperature  without  an  exterm.1  operational  pressure  drop 
as  shown  in  Fig.  B.2a.  The  second  application  is  one  in  which  the  rotameter  is 
operated  at  some  displaced  pressure  and  temperature  with  an  external  operating 
pressure  drop  as  shown  in  Fig.  B.2b.  For  convenience,  the  first  application 
will  be  considered  to  be  the  conditions  for  calibration  of  the  rotameter.  For 
both  applications,  the  bail  position  is  assumed  to  be  fixed,  therefore  K  may  be 
assumed  to  be  constant  if  the  drag  coefficient  is  constant.  If  for  both  applications 
the  fluid  type  remains  constant  and  the  pressure  and  temperature  range  is  small, 
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changes  in  drag  coefficient  may  be  assumed  to  be  small.  For  each  condition, 
equation  (BA)  can  be  used  in  which  the  product  Ky  is  constant  resulting  in  the 
following  two  equations: 

~  -  (' dm!dt)l 

V  (( *-«)*]•» 

(dm/dt)a 

y  i(n-P o)po}> 


Equating  these  two  equations  gives  the  following: 

(dm/dt)j  _  / (ft  -  pi) piY* 
(dm/  dt)0  \(ft-po)p*/ 


If  the  float  is  chosen  such  that  the  ball  density  is  much  greater  than  the  fluid 
density  and  using  dm/di  —  pAV  =  pQ  =  constant,  equation  (B.5)  is  further 
reduced  to  : 


P\Q\ 

P0Q0 


(B.S) 


where 


=  PlQ  2  =  P\Q\ 


Now,  including  the  ideal  gas  lav/  where  p  -  P/(RT),  equation  (B. 6)  can  be 
put  in  terms  of  pressure  and  temperature.  Since  the  drag  coefficient  was  assumed 
to  be  constant  due  to  a  constant  fluid  type,  the  gas  constant  is  also  assumed  to 
be  constant.  Making  these  substitutions  and  solving  for  Q\ ,  equation  (B. 6) 
becomes: 


Q 1 


(B.  7) 
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If  the  temperature,  T\  =  7j  for  adiabatic  conditions,  the  final  correction 
equation  becomes,  for  adie.batic  and  frictionless  flow: 

.5 


-0( 


Po\  (PiTiV 

PcTj 


or  in  terms  of  6P: 


0>.0.((l±i^lWi)J  (B,, 

Generally,  the  calibration  conditions  are  approximate  to  the  operational  con¬ 
ditions.  If  this  is  so,  the  final  equation  will  further  reduce  to: 

Qi  =  Qo  (l  +  (B. 9) 

The  final  equations  (B. 8)  and  (B.9)  can  now  be  used  as  a  means  of  pressure 
correction  the  flow  through  a  rotameter  either  form  a  calibration  curve  or  when 
a  calibration  curve  is  determined  for  which  there  will  usually  be  a  pressure  drop 
associated  with  the  calibration  device.  In  all  instances,  the  calibration  curve 
must  be  determined  and  corrected  for  a  pressure  drop  equal  to  sero.  It  must  be 
noted  that  there  is  a  limitation  on  this  form  the  correction  equation.  The  flow 
from  the  rotameter  to  point  in  which  the  volumetric  flow  rate  must  be  determined 
was  assumed  to  be  adiabatic  and  frictionless.  For  cases  in  which  the  gas  flow 
coming  out  of  the  rotameter  becomes  critical  or  has  heat  added,  there  may  be 
significant  error  to  be  considered. 
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iii)  Rotameter  calibration  and  error  analysis 

Calibration  of  rotameters  is  performed  using  a  calibrated  standard  as  a  ref¬ 
erence.  An  NBS  traceable  calibration  device  is  used  such  as  a  bubble  flow  meter, 
a  dry  or  wet  test  meter,  laminar  flow  element,  etc.  In  addition  to  the  calibrated 
flow  measuring  reference,  a  pressure  reference  to  determine  the  calibration  setup 
pressure  drop  and  additional  calibration  instrumentation  such  as  a  barometer, 
and  thermometer  for  determining  the  ambient  air  pressure  and  temperature  are 
required. 

Each  of  the  above  mentioned  instrumentation  devices  have  some  error  as¬ 
sociated  with  their  used.  This  error  must  be  considered  during  th  calibration 
procedure  to  determine  the  overall  calibration  error  due  to  the  instrumentation. 
Instrument  precision,  accuracy,  and  resolution  must  all  be  considered  when  eval¬ 
uating  instrument  error.  An  additional  error  due  to  the  ability  of  the  calibrator 
of  operator  to  set  the  rotameter  float  at  any  given  setting  consistently  must  also 
be  considered.  This  repeatability  error  must  be  added  to  the  instrument  error  in 
order  to  determine  the  overall  rotameter  calibration  error. 

If  a  function,  ♦,  is  considered  in  which  $  =  ♦(<*,/?, 7),  the  total  error  for  ♦ 
can  be  determined  by  the  following  relation,  (Doebelin,  1983): 

*-i(S9a«+|(S?)a*+i(S?H  <bio> 

where  A  (  )  is  the  error  associated  with  the  measurement  of  the  appropriate 
independent  variable. 
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For  the  calibration  of  a  rotameter,  equation  (5. 9)  can  be  used  with  the 
appropriate  assumptions.  Using  equation  (5.10),  the  error  in  the  calibrated  flow 
rate  can  be  calculated  if  the  instrument  error  in  the  independent  variables  is 
known.  Thus  the  corrected  calibration  flow  rate,  Q„y  for  a  rotameter  can  be 
found  by  the  following  relation  derived  from  equation  (5.9): 


Qo  ~ 


Ql 

(1  +  6P/  5a<m)'5 


where: 


and: 


Ql  =  indicated  calibrated  volumetric  flow  rate  at  an  operational  6P 


Qo  =  corrected  calibrated  volumetric  flow  rate  at  6P  =  0 


the  instrument  error  now  becomes: 

Qo  ! 2AQi  APatm6P _  *P(6P)  \ 

2  V  Qi  (If  6P/Palm)P?lm  (l  +  6P/P.lm)P.tm)  [  •  ] 


_  .  — ,  AV,  V,A< 

5<Jj  =  A<?J  =: 


(5.12) 
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No  /  the  calibrated  volumetric  flow  rate  can  be  determined  with  the  appro¬ 
priate  instrument  error  where  we  have  Qa  ±  Eqq. 

When  the  rotameter  is  in  use,  and  with  the  appropriate  assumptions,  equa¬ 
tion  (15.9)  can  again  be  used  to  determine  the  unknown  Q\  from  the  operational 
pressure  drop,  SP,  and  calibration  volumetric  flow  rate,  Q0.  Now,  the  corrected 
volumetric  flow  rate  at  an  operational  SP  can  be  calculated.  The  method  as 
described  above  can  be  used  to  determine  the  operational  instrument  error  of 
the  rotameter. 

Repeatability  error  can  be  determined  simply  as  the  maximum  deviation  in 
volumetric  flow  rate  when  multiple  repeated  float  settings  are  made.  This  error 
along  with  the  maximum  deviation  due  to  the  instrumentation  error  can  be  added 
to  determine  the  maximum  calibration  error  over  the  range  calibrated.  This  total 
error  can  be  stated  in  terms  of  percent  of  the  range  calibrated  or  in  volumetric 
flow  rate  units  only. 

Actual  use  of  the  calibrated  data  usually  involves  fitting  and  algebraic  func¬ 
tion  to  this  data.  In  this  form,  calibration  information  can  easily  be  interpolated 
between  the  actual  calibration  points.  However,  there  is  also  some  error  associ¬ 
ated  with  the  algebraic  function  when  compared  to  the  actual  calibration  data. 
Goodness  of  fit  tests  and  or  control  charts  with  the  appropriate  confidence  in¬ 
tervals  are  then  requi  -ed  to  determine  the  extent  to  which  the  information  being 
obtained  from  the  calibrated  curve  function  is  correct  or  if  the  calibration  of  a 
specific  instrument  is  remaining  valid.  Generally,  for  a  rotameter,  a  piecewise 
linear  fit  is  sufficiently  accurate  (assuming  there  are  enough  calibration  points) 
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and  within  the  total  rotameter  error  such  that  the  error  due  to  the  curve  fit  can 
be  neglected. 

The  calibration  curve  for  the  rotameter  utilized  in  the  experimented  apparatus 
is  shown  in  Fig.  B.3.  A  dry  test  meter  was  used  as  the  reference  calibration 
standard  with  a  full  scale  error  of  ±0.1  lt/min  over  the  calibrated  range  of  3.0 
lt/min  to  19.0  lt/min. 


B)  Thermocouple  calibration  and  use 
i)  Theory  and  use 

Thermocouples  are  a  first  order  response  temperature  measurement  device 
for  which  dissimilar  metals  generate  a  voltage  (Seebeck  voltage)  which  is  almost 
linear  for  small  changes  in  temperature.  To  measure  temperature  accurately, 
a  reference  junction  must  be  used  in  conjunction  with  the  thermocouple.  This 
reference  junction  generally  consists  of  an  additional  thermocouple  maintained 
at  a  constant  temperature,  usually  that  of  the  melting  point  of  ice  commonly 
referred  to  as  an  ice  point  reference. 

Thermocouples  generally  have  an  error  of  ±1.0,'C  to  ±2.0°C  depending  upon 
the  specific  type  of  thermocouple  used.  If  a  thermocouple  is  individually  cali¬ 
brated  this  error  may  be  reduced  to  ±0.1°C  providing  the  voltage  measurement 
sensitivity  can  be  obtained.  Exposed  junction  iron-constantan  type  thermo¬ 
couples  were  used  lot  these  investigations  requiring  that  the  voltage  measure- 
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ment  sensitivity  be  5.1  /xv  for  a  ±0.1°C  resolution  in  temperature  measurement. 
Switching  between  thermocouples  was  performed  using  an  isothermal  thermocou¬ 
ple  switch  which  was  insulated  to  minimize  stray  emf.  The  thermocouple  output 
was  amplified  with  a  gain  of  100  and  the  resulting  voltage  monitored  using  a 
digital  voltmeter.  The  ice  point  reference  utilized  either  an  ice  bath  or  electronic 
ice  point  for  which  error  induced  due  to  ambient  temperature  variation  was  neg¬ 
ligible.  The  temperature  measurement  setup  was  subsequently  maintained  and 
calibrated  as  a  whole  to  minimize  the  additive  errors  for  each  of  the  components 
such  as  wires,  amplifier,  ice  point  reference,  etc.  An  additional  reference  temper¬ 
ature  to  zero  the  amplifier  was  taken  at  the  ice  point  reference.  This  temperature 
was  periodically  monitored  to  verify  the  existence  of  electronic  drift  and  reset  the 
amplifier  zero  prior  to  taking  temperature  measurements. 

ii)  Thermocouple  calibration 

Thermocouple  calibrations  were  performed  using  an  ERTCO  PG-108  (Ever 
Ready  Thermometer  Co.,  New  York,  N.Y.)  NBS  traceable  (NBS  Identification 
No.  311474)  mercury  thermometer  with  a  resolution  of  ±0.05°C.  Both  the  ther¬ 
mocouple  and  the  calibration  reference  were  mounted  together  and  placed  in  a 
Lauda  circulating  water  bath  and  the  temperature  varied  over  the  range  of  antici¬ 
pated  temperature  measurement  (24.0°C  -  32.0°C).  This  procedure  was  repeated 
and  subsequently  checked  periodically  to  verify  calibration  over  the  course  of  the 
experiment. 

Thermocouple  output  voltage  was  found  to  be  linear  over  the  range  of  tem¬ 
peratures  calibrated  with  a  static  sensitivity,  K*,  of  approximately  5.1  pv/0.1°C. 
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A  typical  calibration  curve  is  shown  in  Fig.  B.4  with  a  99%  confidence  band 
for  the  linear  fit  of  the  calibration  data. 


iii)  Thermocouple  frequency  response 


Exposed  junction  thermocouples  with  a  junction  bead  diameter  ranging  from 
0.025  cm  to  0.076  cm  for  a  time  constant  of  1.35  to  1.73  seconds  respectively 
were  used  for  these  investigations.  Measurement  of  the  transient  temperature 
fields  required  that  information  with  regard  to  the  frequence  response  of  the 
thermocouples  be  determined  to  evaluate  the  attenuation  and  phase  shift  of  the 
measured  transient  temperature  variation.  The  following  relation  for  a  first  order 
instrument  response  shows  the  amplitude  ratio,  AR,  of  the  input  and  output 
signal  or  fraction  of  attenuation  and  the  phase  angle  for  a  given  instrument  time 
constant,  r,  and  frequency,  u>,  of  the  signal  to  be  measured,  (Doebelin,  1983): 


AR  =  —  = 


K' 


L  tan_1(-u;r) 


q%  (w3rJ  +  1) 


.5 


(*■13) 


where: 

AR&  =  20  log  AR 

Figure  B.5  shows  a  Bode  plot  of  a  first  order  instrument  while  Table  B.l 
shows  the  amplitude  ratio  and  phase  shift  for  a  time  constant,  r— 1.7s. 

C)  Hot  wire/hot  film  anemometer  calibration  and  use 


i)  Theory  and  use 


Hot  wire/hot  film  anemometers  consist  of  heating  a  wire  or  film  sensor  by 
passing  a  current  through  a  wire  or  film,  (Doebelin,  1983).  As  a  fluid  passes  across 
the  wire  or  film,  the  amount  of  heat  convected  and  the  resulting  drop  in  sensor 
temperature  for  a  constant  current  or  change  in  the  current  required  to  maintain  a 
constant  sensor  temperature  is  a  direct  measure  of  fluid  velocity.  Typical  hot  wire 
sensors  ate  constructed  of  materials  which  are  certified  to  maintain  consistent 
properties,  the  most  commonly  used  wire  being  a  platinum  coated  tungsten  wire 
with  a  high  coefficient  of  resistance  («0.004  ohms/°C).  Typical  hot  film  sensors 
consist  of  a  fused  quartz  substrate  with  a  high  purity  platinum  film  bonded  to  the 
surface.  Platinum  provides  the  most  stable,  anticorrosive  film  material  available 
with  the  proper  resistance  characteristics  for  thermal  anemometry.  Film  sensors 
are  also  coated  with  either  alumina  for  gas  flows  or  quartz  for  liquids. 

For  equilibrium  conditions,  the  following  relation  relating  the  current  through 
the  sensor,  7,  sensor  temperature,  Tw *,«,  and  sensor  resistance,  72«tr«  to  fluid 
velocity,  V  where: 

I2Rwir«  =  hA  {Twi„  -  Tf)  (2J.14) 

where  the  film  coefficient  has  the  general  form  (King’s  law): 

h  =  C0  +  CiW 


For  these  investigations,  velocity  data  was  obtained  using  a  TSI  IFA-100  hot 
wire  anemometer  using  a  platinum  film  sensor  model  #1211-10  with  a  manufac¬ 
turer  rated  typical  upper  frequency  response  in  sir  at  100  m/s  is  300  kHz.  For 
this  anemometer,  the  sensor  current  is  constant  while  changes  in  the  sensor  tem¬ 
perature  are  measured  by  changes  in  the  wire  resistance.  A  bridge  circuit  is  used 
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to  indicate  a  change  in  sensor  resistance  by  means  of  a  change  in  bridge  voltage 
which  is  directly  related  to  the  fluid  velocity.  The  platinum  film  sensor  that  was 
used  resolves  the  velocity  in  two  dimensions  in  a  cross  flow  configuration.  Velocity 
components  parallel  to  the  sensor  were  not  resolved. 

ii)  Hot  wire  anemometer  calibration 

To  effectively  use  the  hot  wire  anemometer,  the  magnitude  of  the  buoyant- 
induced  velocity,  due  to  the  heating  of  the  sensor,  in  relation  to  the  velocity  to  be 
measured  was  of  interest.  Uncertainty  errors  can  be  determined  depending  upon 
the  orientation  of  the  fluid  velocity  to  that  of  the  induced  buoyant  velocity  as 
described  by  TSI  (TSI  TB  14).  This  is  evident  when  the  stream  velocity  is  opposed 
to  the  sensor-induced  buoyant  velocity  and  of  the  same  order  of  magnitude. 

Calibration  of  the  hot  wire  anemometer  was  conducted  using  a  laminar  cali¬ 
bration  cell  fabricated  to  provide  a  uniform  velocity  profile  as  shown  in  Fig.  B.6. 
Volumetric  flow  rate  was  controlled  and  monitored  using  a  calibrated  rotameter  as 
previously  discussed.  Two  sets  of  calibration  runs  are  shown  (I  and  II)  in  Fig.  B.7 
with  the  calibration  velocity  or  flow  either  being  opposed  (flow  down)  or  aligned 
(flow  up)  with  the  sensor-induced  buoyant  velocity.  This  resulted  in  a  pronounced 
decrease  in  the  output  voltage  for  the  calibration  velocity  being  opposed  to  the 
sensor-induced  velocity  when  compared  to  a  calibration  velocity  being  aligned  with 
the  sensor-induced  velocity.  For  the  specific  sensor  that  was  calibrated  here,  the 
buoyant  induced  velocity  ranged  from  0.010  m/s  to  0.020  m/s  with  a  sensor  sensi¬ 
tivity  ranging  from  0.50  x  10*  -  1.0  x  10J  mv/(m/s)  for  velocities  from  0.020  m/s 


to  0.030  m/s. 
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Discrepancies  between  runs  I  and  II  appear  to  be  related  to  changes  in  the  fluid 
temperature  altering  the  sensor  output  as  would  be  related  to  equation  (i?.14). 
With  a  change  in  fluid  temperature,  the  sensor  resistance  or  temperature  is  shifted 
at  equilibrium  resulting  in  a  shift  in  the  sensor  output  for  the  same  velocity. 


Table  B.l 

THERMOCOUPLE  FREQUENCY  RESPONSE 
AMPLITUDE  RATIO  AND  PHASE  SHIFT  (r=1.7s) 


129 


rad/sec) 

AR 

0.10 

0.986 

0.20 

0.947 

0.30 

0.891 

0.40 

0.827 

0.5C 

0.762 

0.60 

0.700 

0.70 

0.643 

0.80 

0.592 

C.90 

0.547 

1.00 

0.507 

phase  angle  (deg) 

-9.65 

-18.8 

-27.0 

-34.2 

-40.4 

-45.6 

-50.0 

-53.7 

-56.8 

-59.5 
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Fig.  B.l  Free  body  diagram  of  rotameter  float 
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PotQotP* 

6P= 0 


Fig.  B.2  Rotameter  application  Khematic  (a)  operational  preaaore  drop, 
iP  m  0  (b)  operational  preaaure  drop,  tP  /  0 
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Fig.  B.J  rotameter  inlet  air  calibration  curve  corrected  for  SP  — 


operation  with  a  99  percent  confidence  band  for  the  linear  fit 


Fig.  B.5  Thermocouple  first  order  inttrument  Bode  plot 
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Fig.  B.fl  Hot  wire  anemometer  calibration  cell 
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Fig.  B.7  Hot  wire  anemometer  calibration  curves  for  aligned  (flow  up) 
id  opposed  (flow  down)  calibration  velocities 


Appendix  C:  NUMERICAL  FORMULATION  (PHOENICS) 


i)  Ql-File  (two-dimensional  formulation) 

This  is  a  description  of  the  Phoenics  Ql  input  file.  Each  group  is  explained 
as  appropriate  to  convey  the  problem  approach  specific  to  the  Phoenics  software 
package  utilized  for  the  two-dimensional  investigations. 

Group  1:  Run  title  and  preliminaries 

-Problem  title  and  description  with  variable  assignment. 

Group  2:  Time  dependence  and  related  parameters 
-Set  problem  to  transient. 

-Use  method  of  pairs  for  time  step  specification. 

Group  3:  X-direction  grid  specification 
-Set  rectangular  coordinate  system. 

•Default  to  a  one  meter  deep  chamber. 
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Group  4:  Y-dircction  grid  specification 

-Use  method  of  pairs  for  y-direction  grid  layout. 

Group  5:  Z- direction  grid  specification 

-Use  method  of  pairs  for  z-direction  grid  layout. 

Group  7:  Variables  solved,  stored  and  named 
-Single  phase  fluid  specified. 

-Solve  for  and  store  dependent  variables  Pi,  Vl,  Wl,  Hi. 

-Use  implicit  slab-by-slab  method  of  solution  for  velocity  and  enthalpy. 
•Use  implicit  whole-field  method  of  solution  for  pressure. 

•  Use  harmonic  averaging  of  the  exchange  coefficient. 

Group  8:  Terms  in  the  differential  equations 

-Built  in  source  terms  active  for  velocity  components  (pressure  term 


included  in  momentum  equations). 


139 


-Built  in  source  terms  inactive  for  enthalpy  (viscous  dissipation  and 


pressure  work  terms  not  included  in  formulation). 


-Convection,  diffusion,  and  transient  terms  active  for  dependent  variables. 


-Dependent  variables  belong  to  first  phase  in  a  single  phase  formulation. 


Group  9:  Fluid  properties 


-Kinematic  viscosity,  density,  Prandtl  number  specified  for  fluid. 


Group  11:  Variable  initialization 

-Choice  of  restarting  [RESTRT(  ALL))  from  a  previous  solution  or  starting 
from  an  initial  field  specified  using  the  FI1NITQ  command. 


Group  13:  Boundary  conditions  and  special  sources 
-Chamber  inlet  boundary 


Specifies  a  uniform  velocity  and  temperature  on  the  upper  moat  cell  surface 
across  the  width  and  depth  of  the  chamber 
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-Chamber  outlet  boundary 

Specifies  a  reference  pressure  of  zero  across  the  outlet  lower  most  cell  surface. 
This  resembles  the  conditions  similar  to  a  porous  boundary  outlet  in  which  the 
exhaust  side  is  at  a  uniform  and  constant  pressure.  Specifies  an  outlet  enthalpy 
to  be  the  same  as  that  of  the  upstream  neighbor  cells  exhibiting  a  local  parabolic 
behavior  at  the  outlet  boundary. 

-Chamber  vertical  end  wall  boundary 

Specifies  a  constant  temperature  wall  with  no-slip  condition  for  velocity. 

-Chamber  pan  boundary 

Specifies  an  infinitesimally  thick  pan  represented  by  setting  the  velocity  com¬ 
ponents  to  zero  along  the  lower  cell  walls. 

-Chamber  heat  source  boundary 

Heat  sources  are  represented  by  specifying  a  zero  velocity  at  the  cell  location 
within  the  heat  source  area  and  setting  the  temperature  in  the  blocked  region  to 
be  constant. 

-Buoyancy  boundary  using  the  Boussinesq  approximation 

The  Boussinesq  approximation  using  a  ground  subroutine  is  used  to  include 
the  buoyancy  source  term  in  the  vertical  direction.  Note  coordinate  reversal  to 
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be  consistent  with  the  inlet  and  outlet  as  specified  in  Phoenics. 


Group  15:  Termination  of  sweeps 


-Sweep  number  and  residual  limits  are  specified 


Group  17:  Under- relaxation 

-Under-relaxation  is  specified  using  FALSDT  equal  to  1.0.  The  time  step  was 
used  to  control  the  amount  of  relaxation  required  for  specific  solutions. 
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TALK«TSRUN<  1.  1> iVDU-TTY 

GROUP  1.  Run  identifier  and  other  orel  imintne* 

TEXT  (2-D  CHAMBER  FLOW  MODEL  W/  HEAT  SOURCES) 

This  is  a  two-dimensional  model  representing  a  large  rectangular 
chambers  witn  baffles  and  constant  temperature  heat  sources. 

The  dimensions  are  .  £S6m  vertical  and  .813m  horizontal. 


This  is  the  symmetric  transient  model. 

.  VOLUMETRIC  FLOW  (CFM)  - QIN 

.  ....  INLET  VELOCITY  (M/SEC)  - WIN 

.  INLET  TEMPERATURE  (DEG  C)  - TIN 

.  WALL  TEMPERATURE  (DEG  C)  - TWALL 

.  INLET  PRESSURE  (PSIA)  - PIN 

.  INLET  DENSITY  (KG/CU  M) - DENSIN 

_ _  INLET  VISCOSITY  (N  SEC/SO  M>  - VISCIN 

.  INLET  KINEMATIC  VISCOSITY  (SO  M/SEC)  - ENULIN 

.....  COEFFICIENT  OF  EXPANSION  ( 1 /DEG  K)  - BETA 

.  CHAMBER  WIDTH  (METER)  - L 

.  CHAMBER  HEIGHT  (METER)  - H 

.....  NUMBER  OF  TIME  STEPS - NT 

.  TIME  STEP  (SEC) - DT 

.  START  TIME  (SEC)  - TSThT 

.  STOP  TIME  (SEC)  - TSTOP 


REAL.  (QIN,  WIN,  tin,  TWALL,  pin,  DENSIN,  VISCIN,  ENUt  IN,  BETA) 

R£Ai_  (T3TRT,  TSTOP,  DT.  L,  H) 

INTEGER  (NT,  DV1,  DY2.  DY3,  DY4,  DY3,  DY6,  DY7,  DY8,  DYD,  DYIO,  DY1 1 ,  DY12, DY1 3> 
INTEGER ( DZ 1 , DZ2, DZ3, DZ4, DZ3, DZ6, DZ7, A, B, C, D> 

••for  time  steo  C < TSTOP-TSTRT ) -DT*NT3 *• 

DT-. 03 
NT-400 
TSTRT-O. 0 
TST0P-20. O 
PIN-14. 4 
gin-5. a 
TIN-27. 3 
TIN-23.0 
TIN-21.0 
TIN-23. 0 
TIN-24.0 
TWALL— 24. 0 
W I N— Q In/ 1 • 41 14E3 
BETA-l.  0/ <273.  2*TIN) 

DENS IN-PIN* 43. 24/ (492*1. 8«TIN> 

VISCIN-1. 4S8E-6*<273. 2*TIN> •*!. 3/ <383.  6*TIN) 

ENULIN-VISCIN/DENSIN 
L— . 813 

OY 1 -6 1 DY2- 1 | DY3- 1 3  <  0Y4- 1 0 1 DY5- 1 3 1 D Y6- 1 ( D Y7- l 2 
DV8-1 | DY9-1 3 l DY 10- lOjDYl 1-13»DY12-1 I0Y13-6 

ny-Dy l ♦DY2*0Y3*Dv4*Dy3*DY6*0Y7*DY#*DY9*DY 10*DY ll*DY12*DY13 

A-DY1*DY2*DY3*0Y4*DY5*0YB 

B-NY-DY13-DY12-DY1 1 -DY10-DY9-DY8*! 

C-NY-DY 13-DY12-DYI 1*1 
h— , £86 

DZl-38|DZ2-8|DZ3-£|DZ4-33|DZ3-8|DZB-2|027-7 
NZ-DZ  1*0Z2*0Z3*DZ4*DZ.3*DZ6*DZ7 
D-nZ-0Z7-0Z6-0Z3*1 
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GROUP  2.  Tim*  deoanoenc*  and  related  parameters 
STEADY-F 
TFIRST-TSTRT 
Tl_AST=*TSTOP 
LSTEP-NT 
TFRAC  <  1  )  —NT 

TFRAC (2) -DT/ ( TSTQP— TSTRT ) 

GROUP  3.  X-flir«ctlon  grid 
CARTES-T 

GROUP  4.  Y-d  l  r-*Ct  1  on  grid 
YVLAST  =i_ 

YF  RAC ( 1 ) — DY1 ; YF  RAC (2) *6. 0E-2/DY1 
YFRAC (3) »DY2i YFRAC (4) -1. OE-2/DY2 
YFRAC(5> »DY3? YFRAC (6) ». 13/DY3 
YFRAC  ( 7  >  «DY4 ;  YFRAC  ( 8 )  *.  10/DY4 
YFRAC (9) »DY5 ( YFRAC < 10) «. 13/DY5 
YFRAC (11) *0Y6i YFRAC  t 12) -1. OE-2/OY6 
YFRAC( 13) -DY7» YFRAC < 14>«. 12/DY7 
YFRAC  ( 15) «DY8l YFRAC  ( 16) -1. 0E-2/DY8 
YFRAC < 17) -DY9 1 YFRAC  < 18) 13/DY9 
YFRAC (19) -DY 10 » YFRAC (20) «. 10/DY10 
YFRAC (21 ) -DY1 1 ; YFRAC (22) «. 13/DY1 1 
YFRAC (23) -DY12; YFRAC (24) -t . 0E-2/DY12 
YFRAC (25) -DY13| YFRAC (28) -6. 0E-2/DY13 
GROUP  5.  Z-dir*ction  grid 
ZWLAST-h 

ZF RAC  < 1 ) — DZ 1 i ZFRAC (2) *. 38/DZ1 
ZFRAC (3) *DZ2 l ZFRAC  <  4 ) ■S. 0E-2/DZ2 
ZFRAC (5) “DZ3| ZFRAC (6) »2. 0E-2/DZ3 
ZFRAC (7) -DZ4, ZFRAC (8) 3S/DZ4 
ZFRAC (9) «DZ5i ZFRAC < 10) -8. 0E-2/DZ5 
ZFRAC (11) «DZ6«  ZFRAC ( 12) -2. 0C-2/DZ6 
ZFRAC ( 13) -DZ7, ZFRAC ( 14) -7. 06-2/DZ7 

GROUP  8.  Body-fitted  coordinate*  or  grid  dittortion 
GROUP  7.  Varuolti  tolvtd,  itortd  ana  named 
0NEPhS»T 

SOLUTN  <  PI , Y, Y, Y, N, N, Y) 

S0UUTN(W1  , Y, Y, N, N, N, Y) 

S0LUTN(V1,  Y,  Y,  N,  N,  N,  Y> 

SOLUTN (HI , Y, Y, N, N, N, Y) 

GROUP  8.  Term*  (in  differential  Kuationt  and  device*) 

TERMS (VI , Y, Y, Y, Y, Y, N> 

TERMS  <U1  ,  Y,  Y,  Y,  Y,  Y,  N) 

TERMS (Ml , N, Y, Y, Y, Y, N) 

GROUP  9.  Prooerti**  of  tn*  medium 
ENUL"ENULIN 
RH01-DENSIN 
PRNDTL (HI ) -. 7 

GROUP  10.  1 nt*> — Bni»*-tr#n*f*r  oroce****  and  oroeertie* 

GROUP  11.  1 n 1 1 i a  1 i rat i on  of  variable* 

FIINlT(Hl)-( T in* T WAUL ) /2 
F11N1TIWH  -O.  0 
FI  I N 1 T (VI )-U. 0 
F I  I  NIT (PI > »0. 0 
REST  RT (ALL) 


144 


GROUP  12.  Convection  and  icn  JO 

GROUP  13.  Boundary  CC-naitlOfi*  ar.c  latcui  source* 

••inlet  3our.Oiry  ecnflltlOft** 

PATCH ( INLET, LOW. 1, 1, 1, NY. 1, I,  l.*,7> 

ccvGl.  <  :nlet,  pi,  -  i  *flu.  rho i »win> 

COVAl  (  I NLCt .Hi,  QNL  YHS,  W  I  N  ) 

CCVflt < INlET, - : . ONLYHS, T IN' 

•♦outlet  SOunoary  coriflitiori** 

patch (Outlet,  high, i , i, i, ny, nz. nz. :.  n~> 

COVAl ( OUTLET ,P1,FIXP,  0.  O) 

COVAl  (CuTwET,  HI,  QN-.yWS,  SAHE) 

**••*11  Boundary  condition*** 

PATCH  (  mAl  _  1 ,  ShAll,  t,  1,  1,  1,  1 .  NZ, 

COVAl  (hALLI,  hi,  1.  A/PRNDTL  <H1  i ,  TWA..,.) 

COVAl (wAllI, HI, 1. 0, 0. 0) 

COVAL  ( MALL  1 ,  V 1 ,  1.0,  0.0) 

PATCH  (HALl2,  NWALL,  1,  1  ,  NY ,  NY ,  1,NZ,  1 .  NT  > 

COVAL  (  WALL2,  Hi,  1. 0/PRHDTL(Hl  >  ,  TWAUu  > 

COVAl (WALL 2, HI  ,  1 . O, 0.  0> 

COVAL  <WAll2,  VI,  1. 0,  0.  O) 

**o an  rtttnetion  represented  By  setting  v*.  :eity»i.i*» 

PATCH < PAN 1, HIGH, 1, l, DV1*1, A, DZl*DZ2*DZ 3.  DZ l*DZi*DZ3, l, NT> 
coval <pani,wi,fixval,o.o> 

COVAl*'  .  VI,  FI  XVAl,  0.0) 

PATCH  hIGH,  1,  l,B,NV-0Vl3,D2l*DZ2OZ3.DZ:*0Z2*DZ3.  I, NT) 

COVAl  ..,  HI,  FI  XVAl,  0.  0) 

COVAl  <  vAN2, VI , F I  X VAL,  0.  O) 

PATCH <PAN3, HIGH,  1,  1,0V 1*1, A. NZ-DZ7.  LZ-0Z7.  1,N") 

COVAl (PAN3, HI, FIXVAL,  0.  0> 

COVAl  (  PAN 3,  VI  ,  F  I  XVAL,  0.  0) 

PATCH (P4Na, HIGH,  1,  1,B, NY-0V13,  NZ-DZ7,  NZ-0Z7,  :,  NT) 

COVAL  (PAN-*,  Ml  .  F  I  XVAl, O. O) 

COVAL  <  PAN4,  VI  ,  F  I  XVAL,  C.  0) 

•  •r>eat  lO'jTct  reoresentat  ion** 

••Heater*  1*« 

PATCH (HTR1. CElL,  1,1, DV1*Dv2*i,Dyi*Dy2*Dy3,  DZ 1 ♦ 1 ,  DZ I *DZ£,  1,  NT) 

COVAl  <hTB1 , HI , F I XVAL, O. 0> 

COVAl (hTRI, VI, FI XVAL, 0.  O) 

COVAL (HTRl , HI . F I XVAL, 32. 0) 

••neater  £•• 

PATCH (HT»a, CELL, 1,  1, DV1*0Y2*DY3*DV4*!,  a- DVB,  DZ I  *  1 .  DZ 1 ♦DZ2,  1, NT) 

COVAl  (HT  R«r,  HI  ,  F  I  XVAL,  V.  O) 

COVAl  (HTR2,  VI,  FI  XVAL,  0.0) 

COVAL  <HT»£, hi, F I XVAL,  32.  0) 

••neater  3** 

oorr.MiHToi,  cell.  I.  I,  B-rDY#,  NV-Dv  13-D  V  U  -  Dv  j  1  -Dvi<>,  DZ  1*1.  DZ1*DZJ',  1 .  NT  > 
COVAl  IMTOJ,  HI  ,  F l XVAL,  O. O) 

COVAL (HTR3, VI , F  J  XVAl, 0. 0> 

COVAL  <hTR3, HI , F I XVAl, 32.  0) 

••neater  »•• 

patch ( MTU*. CEwL,  1, ;,C, NV-DV 13-Dvii.  oz 1*1.  DZ.*D Z2.  t.sT) 

COVAL  <wTB4, hi , F I XVAL, O.  O) 

COVAL  CHT*4, VI, F I XVAL.O. 0) 

COVAL  (HT*4,  hi  ,F  1  XVAL,  32.  0) 

••neater  3*e 


< 
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PATCH <mtP5. CEll,  1,1, DYl ♦DV2-M  , DYl  ♦Dy2«-DY3, 0, NZ-DZ7-DZ6,  1 , NT) 

CCVAl  (HTR3,  wl.FI  xvAl,  0.  0) 

COvAl (hTBS. VI.FI xvAl, 0. 0) 

COvAl  ih-B3,  Hi  ,  F  I  x  v Al ,  32.  0> 

*»n**t#r 

PATCh(ht»6.C£-_,  .  1 ,  CY  :  ♦DY2»DY3*DY4*  l ,  A-DY6,  D,  NZ-DZ  7~0Z6,  1,nT) 

C  C-V«h!_  <  mt  q£> ,  *  V&L  «  C  .  ‘ 

~  2  i  ~  .  ✓  1  .  *  \  '  -  .  ■  > 

C  Z  vkL  <  h'1'  ►"  1 .  c  «  *  -  *-«)-  ,  2  w .  0 ) 

**n**t*r  7»* 

PATCh <hTR7. C£ _ _  1,1. 6*Dva. ny-Dyi3-Dy 12-DY1 1-OYlO,  0,  NZ-0Z7-DZ&,  l,  NT) 

COvAl (mTR7,  wl.  F I  X  V Ac , O. O) 

COVAL  <mT«7, VI.  FI  XVAL,  0.  O) 

COvAl  (hTR7.i«1,FI  X  VAl,  32.  0) 

8«* 

PATCH  ;htB8,  CEl._,  1.  1.  C.  NY-DY13-0Y12.  D.  NZ-D27-DZ6,  1.  NT) 

COVAl (hTB®,  Nl.FI xvAl, 0. 0) 

COVAl (HTPa. V1.fi XVAL. O. 0) 

COvAl  <  H-oa. Hi . F I x vAl,  32.  0) 

•  •Ouoyxwcy  SOunfliry  —  E»Ou»»|H««a  *  O  pro*  1  W4t  1  on  *• 

•  ••ourct  t#rw  -  BS3 1  fi*  (  BSG  1  ♦BSG2*h  1  )  or  g»o»t  •  <  T  j  n-r  )  •* 

••(not*  ccorflinxt*  r#v*r**l>  ** 

hsoio»9. ai 

PATCH  (PuOY,  PHASE «*.  1.  1.  1.  NY.  1.  NZ.  1  .  NT  ) 

COVAL (PuCy. Nl.FI xFlu, 3BND3 ) 

«S32»-E>£’A 
FfSGl-BETA  •  -IN 

GBGUP  14.  Dc— r#«M  o rNlur*  for  PABCB-.  T»Ut. 

GBGlP  15.  ’»r«lRltlO«  Of  l«NOI 
L ShEEP-1" 

E-fc 

£  a 

e-t. 

,  e-a 

’trxi'XtlOO  o f  IO"» 

•  1 *<«t  ion 


P£  5BE  F  (PI ) «1  . 

OESOEF (HI ) -1 . 

OES»EF  < VI ) «1 . 

BE  S  BE  b (Ml  )•! , 

g»olp  is. 

&O0VP  1  7 .  Ur.(M 
»E  t_AX  (Ml  ,  CA.  30T.  1 .  '» 

«Eii_AX<Vl.FA^SD',  1  .  •'»> 

OEwAX  (M1.FOL50T,  1.0> 

BEL  a*  « oi . f&v SDT,  1 . 0» 

GOCuP  I*.  Ll*)t«  O"  vI’-IIOlRI  or  incrNMl  to  t'#* 
r,»CvP  I  a.  r«t4  CO— .p  IC«t  rn<3  By  t*l*U  It*  to  3*0tjn*0 
'jOOvA  20.  Pr*J  i-i  »  «ry  print-ou* 
f  C.Hj*F 

G»0l  D  >1.  F  l  #  ;  ij  ar.O  — onitor  COOlTdl 


or  1  nx  *  O',  t 


-M  ;1Pi  T  l  P  1  .  »  . 

S,  S,  N,  *, 

v  » 

„  *  Pi )  *  x  l  .  y  . 

s.  s,  s.  *, 

*  > 

)■  /  - 1»  ’  'Vi,  y, 

S,  N,  N,  y, 

*  i 

Ov’BVT  (h1  .  y. 

n,  s,  s,  y, 

*  > 

o  t 

c  .  ••  r>  ■ '  * 

V  S  l 

nOOMON- 

I  i  »'  >.*•,  >  /  £ 

I  l  "O’hN  l  '  4 
!  PL  *L*l  VwE  E  p 

GB'aJB  f  t. 


F|«1R  C»r(oX*<9ut  BlOt  CC«tXOl 
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NPLt-1 

css: z». 4 

\t'ORIN«NV/30 

nzp»in«*z / 10 

NTPS IN«N~/5 

I  T 

STOP 


ii)  Ql-File  (three-dimensional  formulation) 


This  is  a  description  of  the  Phcenics  Ql  input  file.  Each  group  is  explained 
as  appropriate  to  convey  the  problem  approach  specific  to  the  Phoenics  software 
package  utilized  for  the  three-dimensional  investigations. 

Group  1:  Run  title  and  preliminaries 

-Problem  title  and  description  with  variable  assignment. 

Group  2:  Time  dependence  and  related  parameters 

-Set  problem  to  transient. 

-Use  method  of  pairs  for  time  step  specifi’ation. 

Group3:  X-direetion  grid  specification 

-Set  rectangular  coordinate  system. 

-Unifoun  grid  spacing  io  the  a-dirrction  with  five  cells  for  a  chamber  depth 
of  6.668  cm. 

Group  4:  Y  direction  grid  specification 
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-Use  method  of  pairs  for  y-direction  grid  layout. 

Group  5:  Z-directioa  grid  specification 

-Use  method  of  pairs  for  z-direction  grid  layout. 

Group  7:  Variables  solved,  stored  and  named 
•Single  phase  fluid  specified. 

-Solve  for  and  store  dependent  variables  Pi,  VI,  Wl,  Hi. 

•Use  implicit  slab-by-slab  method  of  solution  for  velocity  and  enthalpy. 

•  Use  implicit  whole-field  method  of  solution  for  pressure. 

•Use  harmonic  averaging  of  the  exchange  coefficient. 

Group  8:  Terms  in  the  differential  equations 

-Built  in  source  terms  active  for  velocity  components  (pressure  term  included 
in  momentum  equations). 

•Built  in  source  terms  :  a  act  ire  for  enthalpy  (viscous  dissipation  and  pressure 
work  terms  not  included  in  formulation). 
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-Convection,  diffusion,  and  transient  terms  active  for  dependent  variables. 


-Dependent  variables  belong  to  first  phase  in  a  single  phase  formulation. 


Group  9:  Fluid  properties 

-Kinematic  viscosity,  density,  Pranatl  number  specified  for  fluid. 


Group  11:  Variable  initialization 

-Choice  of  restarting  [RESTRT(ALL)J  from  a  previous  solution  or  starting 
from  an  initial  field  specified  using  the  FIINITQ  command. 


Group  13:  Boundary  conditions  and  special  sources 
-Chamber  inlet  boundary 

Specifies  a  uniform  velocity  and  temperature  on  the  upper  most  cell  surface 
across  the  width  and  depth  of  the  chamber 

-Chamber  outlet  boundary 

Specifies  a  reference  pressure  of  sero  across  the  outlet  lower  most  cell  surface. 
This  resembles  the  conditions  similar  to  a  porous  boundary  outlet  in  which  the 
exhaust  side  is  at  a  uniform  and  constant  pressure.  Specifies  an  outlet  enthalpy 
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to  b«  the  same  as  that  of  the  upstream  neighbor  cells  exhibiting  a  local  parabolic 
behavior  at  the  outlet  boundary. 

-Chamber  vertical  end  wall  boundary 

Socdfies  a  constant  temperature  wall  with  no-flip  condition  for  velocity. 

•Chamber  vertical  side  wall  boundary 

Specihes  either  a  constant  temperature  wall  or  defaults  to  an  adiabatic  wall 
with  no-slip  condition  for  velocity. 

-Chamber  pan  boundary 

Specifies  an  infinitesimally  thich  pan  represented  by  settii  g  the  velocity  com¬ 
ponents  to  sero  along  the  lower  cell  walls.  Pane  spaa  the  depth  of  the  chamber. 

-Buoyancy  boundary  using  the  Boussinesq  approximatioc 

The  Eoussinesq  approximation  using  a  ground  subroutine  is  used  to  include 
the  buoyancy  source  term  in  the  vertical  direction.  Note  coordinate  reversal  to 
be  consistent  with  the  inlet  and  outlet  as  specified  in  Phoenics. 


Group  15:  Termination  of  sweeps 


-Sweep  number  and  residual  limits  are  specified 
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Group  17:  Under- relaxation 


-Under-relaxation  is  specified  using  FALSDT  equal  to  1.0.  The  time  step  was 
used  to  control  the  amount  of  relaxation  required  for  specific  solutions. 
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TALK-TsRUNt  1.  1);VDU«TTY 

GROUP  1.  run  id»nt»fi*r  and  otn»r  orthrmnarm 
TEXT (3-D  CHAMBER  FLOW  MODEL  W/0  HEAT  SOURCES) 

Tni*  is  a  thr«#-a i m«n*i onal  mod»l  rtarttant  mg  ths  laaorwtory 
iKDtnmtntal  wostl  mtn  Oaf^le*. 

"r-i»  dimensions  are  .636m  vtrtical  and  .SI 3m  nerizontal 
anu  6. 67cm  aeeo. 

T11S  IS  t"e  symmetric  transient  mods!. 

.  VOLUMETRIC  F_3w  (CFM)  - QI\ 

.  Inlet  VELOCITY  (m/SEC)  - - — — - —wiN 

.  INLET  TEMPERATURE  (DEG  C)  - TIN 

.  wAt_L  TEMPERATURE  <0E3  C)  - TWAll 

.  INLET  PRESSURE  (PSIA)  - PIN 

.  INLET  DENSITY  (KG/CU  M) - DENS  IN 

.  INwET  VISCOSITY  (N  SEC/SO  M)  - VISCIN 

.  INLET  KINEMATIC  VISCOSITY  (SO  M/SEC)  - ENULIN 

.  COEFFICIENT  OF  EXPANSION  ( I /DEG  K>  - SETA 

.  SPECIFIC  HEAT  PHASE  1  (AIR)  (J/KG  DEG  K)  —CPI 

.  wall  thermal  conductivity  <w/n  DEG  C>  - kcqnw 

.  wall  thickness  (m> - - - *#alth 

.  NUMBER  OF  TIME  STEPS - NT 

.  T I  ME  STEP  (SEC)  - DT 

.  START  TIME  (SEC)  - TSTRT 

.  STOP  TIME  (SEC)  - TSTOP 

REAL (QIN.  WIN,  TIN, TWALL, TROOM.KCONW, PIN,  DENSIN,  VISCIN,  ENULIN, BETA) 

REAL ( TSTRT. TSTOP, DT, CPI, WALTH) 

INTEGER (NT) 

••for  ti me  steo  C (TSTOP-TSTRT) -DT*NTI 
DT*. 3  ' 

NT* 1 00 

tstrt»o. o 

TSTOP-50. O 
CP  I *100*. 0 
PI N-14. * 

uin»3. a 
’IN-37. 3 
TIN-35.  (» 

’’IN-31.  O 
TIN-33. 0 
*  I N-34. O 
0 

*  ROCM-ce. 73 
WiN-OTM/J.  41  t  *E3 
EETP-1. 0/ (373. 3-tjn) 

DEN3IN-PIN-43. £4/ (433*1. a*TIN) 

v: sc  IN-1 , 43BE-6* (£73. 3- TIN) •*!. 3/(363. 6*TJN> 

ENUL IN-VI  SC  IN/ DENS  IN 
L-.8U 

DY 1 -6 1 DY3-1 i DY3-| Ji DY4- 1 Oi OvS- 1 3 1 DY6- 1 |DY7- 1 2 

DY8-1 i  Dy*J-  i  3 »  Dy  1 0- 1 0 1 Oy i 1  - 1 3 1 Dy 1 3- 1 1 DY 13-4 

NV-DYl *DV3*DY3»DV4*DY3*Oyfc*DY7*OYa*DYt*OVtO*DYl i *DYJ3*0Y13 

A-Ov i *Dy3*0y3*Dy4*0y3*0yS 

B-NY-0Y13-DY13-DYI 1 -OY10-OYS-DVi- 1 
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C»(ny-DV  12-DY  i£-DY  1 

h*.  e,a& 

DZl=*28;DZCr=d:  OZo-d:  s  DZ 4-35  s  DZ 5-8  s  DZ6-d  *  OZ  7-7 
NZ-DZ1  ■OZivDZ20Z4  +  DZ5*DZ£*DZ7 
D-NZ-3Z7-DZb-DZS-l 

GROUP  Z.  time  aeoer o»nct  ar.c  related  oararneter* 
5-EOD  <  =  p 
-f : rst  =  "t  =  ’ 

T_-5’=T3~3P 
.  ■  ’’i 

'  p-^AC  <  1  )  »-nT 

--RfiC(Z-)  »DT/  (  TSTGP-TST  RT ) 


SRCuP  3.  x-cxrect ion  grid 
Cfi-STzS«T 

3RDP-R<x.5.  6.  fctdE-,2.  1.0) 

GRCUP  4.  y-direction  grid 
YVtOST-u 

Y- ROC ( 1 ) «-DYJ ; YPROC  <2> -6. Ofi-S/DYl 
YPROC  <  2) «DY2  « VP ROC (4) -1. 0E-S/DV2 
» - RQC < 5) *DY2 j v=eoC <&) «. 12/0Y3 
vc  ROC (7) -Dv4 :  YPROC < 8 )  «.  10/DY4 
i-  ROC  O >  «P  v5 i vc  ROC (  10)  -.  1 3/DyS 
if  ROC  (  1  1  )  *■  DYfcj  vPROC  < ld> -1.  CE-i/DY8 
*  *  ROC < 1 3) -Dv7( ye  ROC <14  1-.  li/0 Y 7 
Y=R0C  < 15) -JY8» VP  ROC  <  1  fc>  -1. 0E-i/DY8 
v- ROC  < : 7 ) -0v3i vproc (18)-.  12/DV3 
yp  ROC  t  1  3'  Ov  10|  vP p’OC  <i0> -.  10/ Dv  10 
v : ROC  <2 1 ) «DV : : i vP  ROC <ii> -.  13/DY11 
vpROC  <  i  3  >  -  Dv  l  i  i  yp  ROC  <i4>  -1 .  0E-£/DY li 
» -  ROC  <  £5  >  “Dv  i  j  i  vp  ROC  <i6>  -fc.  0£-i/DYi3 
GROUP  *.  1-directior  grid 

Z-L&Sr-H 

Zc  ROC < 1 ) — OZl  i ZPROC  <i> 38/ DZ 1 
Z P °0C  <  3 )  «DZ<  i  ZP  ROC  <  4 )  -8.  Ofe-i/DZi 
Z'-'“OC  15)  -DZ2i  ZPROC  <fe)  -d.  0E-2/DZ3 
ZPROC  <  7) -DZ* : ZPROC (8) -.  35/D Z 4 

:■  -oc  <  ji  »oz*  i  z p -oc  <  iO)  -8.  <>e-i/oz5 

ZPROC  I  1  1  )  -028  i  ZPROC  <  Id)  -i.  O€-i/0Z8 
ZPROC  <  1  3  )  -02  /  l  ZcROC  <  1  4  )  »-7.  >.'£-d/DZ  7 

GROUP  8.  fccdy-fitted  cocsinttf*  or  grid  oittortion 
GRQuo  7,  vart«Dle«  •olvea,  stored  and  named 
Q\£PhS-t 

•CL  J7N<  PI  ,  V,  V.  V,  N.  N.  V) 

CC^UTN <  W1  ,  V,  V ,  N,  N,  N,  V ) 

1  ?.  U  TN  <  VI  ,  V,  V,  N,  N,  N,  V  ) 

CCl  ;tN  <U-.  ,  v,  v,  N,  n.  v> 

,  T  V  I  -  1  .  V  .  V,  N,  S,  v  ) 

i"i»0up  8.  term*  i  in  oi<<f  fntii.  eo'.iaticn*  ana  device*) 
'  E  R  -  3  (  V  l ,  V,  v,  v,  v,  y.N) 

• w- c,  Mjl.v.y.v.v,  v,N) 

'<»*!)  <  wt ,  v,  y,  y ,  V,  V,  N) 

’  £  R«is  <oi,n,  v.  y,  y.  v.s) 

CiROup  3.  o»*r*oer  t  i  e*  of  t*e  -ed  i  <jn 

NUL  *fe  NUL  IN 
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RH01-DENSIN 
PRNDTL  (HI  )  ■.  7 

GROUP  SO.  Inte> — w*»t-tr*n*f«r  orocnm  and  Drooert i e* 
GROUP  11.  initialization  of  variables 
FI  INIT  (HD-TROOM 
FI INIT (W1 ) »0. 0 
FI INIT (PI ) *0. 0 
RESTRT (ALL) 

GROUP  IS.  Convection  and  diffusion  adjustments 
GROUP  13.  boundary  conditions  and  soecial  sources 
♦•inlet  boundary  condition*# 

PATCH  < I NUET , LOW,  1, NX, l, NV, 1, 1* l»N7) 

COVAL ( INUET, P1,FIXFLU, RHQ1*WIN> 

COVAL  ( I NUET, Ml, ONLYNS,  WIN) 

COVAU ( INLET, HI , ONLYMS, TIN) 

•♦out  1 et  boundary  cond 1 1 1 on** 

PATCH (OUTLET, HIGH, 1, NX, 1,  NY, NZ,  NZ,  1,NT> 

COVAL (OUTLET, P1,FIXP,  0.0) 
coval (Outlet, hi , onlyms, same) 

••wall  boundary  conditions** 

PO*CH( WALLS, Shall,  1,NX,  1,  1, 1, NZ.  l.NT) 

COVAL  <  WALL ) , H 1 ,  1 . 0/ PRNDTL  < H 1 ) , TWALL ) 

COVAL  (WALL  1,  Wl,  1.0,  0.0) 

COVAL (WALL1, Ul, 1. 0, 0. C) 

COVAL (WALL  1, VI, 1.0,  0.0) 

PATCH (WALL2, NWALL, 1,NX,NY,NY, 1,NZ, 1,NT) 

COVAL (WALLS, Ml, 1 .  O/ PRNDTL < Ml >, TWALL) 

COVAL ( WALLS, Wl , 1.0,0. 0) 

COVAL (WALLS. Ul,  1.0,0.  0) 

COVAL ( WALLS, V l , 1 . 0,  0.  0) 

PATCH(WALL3, WWALL. 1,  1,  i, NY, 1,NZ, 1,  NT) 

COVAL (WALL3, Ml, 1 . 0/ PRNDTL (Ml )  ,  TRQQM) 

COVAL (WALL3, Wl, 1.0,  C.O) 

COVAL  ( WALL 3,  Ul,  1.0,  0.0) 

COVAL (WALL3, VI,  1.0, 0.0) 

PATCH (WALL*, EWALL, NX, NX, i,NV,  1,NZ, 1,NT> 

COVAL (WALL*, Ml,  1. 0/ PRNDTL (Ml ) , TROOM) 

COVAL (WALL*. Wl, 1.0,0. O) 

COVAL (WALL*. Ul,  1.0,  0.0) 

COVAL (WALL*, VI , 1 . 0, 0. O) 

•*oan  restriction  reoresented  by  setting  velocity"^)** 
PATCH  ( RANI ,  HIGH,  1,  NX,  DY  1*1,  A,  DZ  l*DZS*DZ3,  DZl*DZS*DZ3,  1,  NT) 
COVAL  (  PAN  1 ,  W 1 ,  F  I  X  VAL  ,0.0) 

COVAL ( PAN1,  Ul, FI XVAL,0. 0) 

COVAL (PANl,  VI, FIXVAL, 0. 0) 

PATCH (PANS, HIGH, l, NX, B, NY-DY13, DZI*DZS*DZ3,  DZl*0ZS*DZ3,  1,  NT) 
COVAL ( pans, Wl , F I XVAL, O. O) 

COVAL ( PANS, U1.fi XVAL. 0. 0) 

COVAL  (PANS,  VI,  FI  XVA<.,  0.  0) 

OATCH(PAN3,  HIGH,  1  ,  NX,  DY1*  1,  A,  NZ“DZ7,  NZ“DZ7,  l, NT) 

COVAL ( PAN3.  Wl , F l XVAL, 0. 0) 

CCVAL  (  PANJ,  Ul  ,  F  1  XVAL,  0.  O) 

COVAL (PAN3, V) , FI XVAL,  0.  0) 

0ATCH(PAN*, high, i, NX,  B, NY-0Y13, NZ-DZ7,  NZ-DZ7,  1,  NT) 

COVAL ( pan*, Wl, Ft XVAL, 0.0) 
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COVAL (PAN*, Ul, FIXVAl, 0. 0) 

COVAL (PAN*, VI, FIXVAl,  0.  0) 

»«oouy*ncy  Boundary — (ouitintig  aooromngtton  ** 

•  ♦sourct  term  ■  RSG 1 0* < RSG 1 ♦ RSG2*H1 )  or  g*bet a ( T i n-T)  ** 
*♦ (not*  coorain*t»  reversal)  •* 

RSG 10-9. 31 

PATCH (BUOY, PHASE*, 1, NX, 1. NY, 1, N2, 1 . NT ) 

COVAL (BUOY, ui, ri xflu, grnd3> 

SSG2--BETA 
R33 1  * BE T A  *  TIN 

GR2LP  1*.  Doir.itriim  oressure  for  PARAB*».  TRUE. 

GROUP  13.  termination  of  immoi 
lSWEEP-10 
RESREF ( Pi > ■ 1 . E-6 
RESREF (W1 l-l.E-6 
RESREF (Ul > -t. £-6 
RESREF (VI >-l. £-6 
R£SREe (HI ) *1 . E-6 

GROUP  16.  Termination  of  iteration* 

GROUP  17.  unoer-rel a«*t ion 
RELAX (UI.FAlSDT,  1.0) 

RELAX (UI.FAlSDT, 1.0) 

R£_AX (VI, FAlSDT, 1.0) 

RElAX (HI, FAlSDT, 1.0) 

RELAX (PI, FAlSDT,  1.0) 

GROUP  IB.  Limits  on  vanaole*  or  increments  to  them 
GROUP  19.  Data  communicated  by  satellite  to  GROUND 
GROUP  20.  orelinmary  orint-out 
ECMO-F 

GROUP  21.  field  and  monitor  control 
OUTPUT (Pi , Y, N, N,  N,  Y,  Y> 

OUTPUT  (Ul,  Y,  N,  N,  N,  Y,  Y) 

Output  <ui , y, n, n, n, y, y> 

Output  <  vi,  y,  n,  n,  n,  y,  y> 
output (HI  ,  Y,  N,  N,  N,  Y,  Y) 

GROUP  22.  soot  value  orint-out 
npbhOn-20 
J ymON-30 
I ZHON-23 
!  xr-ON-3 

:p_tl-l sheep 

GROUP  23.  field  orint-out  and  Plot  control 
NPL  T • 1 

c«si * 

NYPRIN«3 
nZPRJN«NZ/ 10 
NXPR jn«3 
NTPR1N-NT/2 
I *PBl«3 
STOP 


Appendix  Dj  DISCUSSIONS  ON  THE  USE  OP  PHOENICS 


Using  "PHOENICS”  for  this  investigation  involved  extensive  computational 
time  due  to  the  physical  size  of  the  problem.  A  100  x  100  grid  was  used  in  the 
solution  of  a  problem  with  the  momentum  and  energy  equations  being  coupled 
through  the  buoyancy  body  force  term.  In  addition,  the  time- dependent  formula¬ 
tion  was  chosen  to  address  time- dependent  behavior  adding  to  the  computational 
time.  The  computational  time  was,  however,  reduced  by  increasing  the  F  array 
and  NFDIM  in  the  ground  from  a  default  of  50000  to  300000.  Increasing  these 
dimensions  effectively  decreased  the  I/O  to  the  disk  and  increased  amount  of 
operation  within  the  core  memory.  No  additional  modifications  to  the  ground 
subroutine  made  for  this  investigation. 

The  buoyancy  body  force  was  incorporated  using  a  ground  subroutine  which 
utilises  the  Boussinesq  approximation.  The  only  difficulty  in  the  specification  of 
this  ground  subroutine  was  the  direction  of  the  body  force  term.  With  "PHOEN- 
ICS",  the  inlet  along  the  s-axis  is  taken  to  be  the  "low”  with  the  outlet  being 
the  "high".  A  coordinate  reversal  was  required  to  maintain  consistency  with 
the  direction  of  the  body  force  term  and  the  direction  of  flow  in  the  s-direction 
specific  to  "PHOENICS". 
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